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ABSTRACT 


The  success  of  interceptors  used  by  the  United  States  ballistic  missile  defense  program  is 
jeopardized  by  the  use  of  hostile  missile  decoy  deployment  and  evasive  maneuvers.  The 
ability  to  discriminate  between  legitimate  threats  and  decoys  is  a  crucial  requirement  for 
interceptor  algorithms.  The  feasibility  of  the  interacting  multiple  model  joint  probabilistic 
data  association  filter  to  effectively  track  a  ballistic  missile  and  detect  decoys  and 
maneuvers  is  the  focus  of  this  thesis.  Model  development  and  data  association  schemes 
are  discussed  along  with  optimized  values  for  selected  parameters. 

Performance  comparisons  of  the  resultant  algorithm  to  a  standard  Kalman  filter 
utilizing  a  nearest  neighbor  discriminator  are  conducted.  Scenarios  include  combinations 
of  missile  maneuver  and  decoy  deployment.  While  the  Kalman  filter  experiences  limited 
success,  the  proposed  filter  tracks  the  missile  in  every  scenario. 
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EXECUTIVE  SUMMARY 


The  success  of  interceptors  used  by  the  United  States  ballistic  missile  defense  program  is 
jeopardized  by  the  use  of  hostile  missile  decoy  deployment  and  evasive  maneuvers.  The 
ability  to  discriminate  between  legitimate  threats  and  decoys  is  a  crucial  requirement  for 
interceptor  algorithms.  The  feasibility  of  the  interacting  multiple  model  joint  probabilistic 
data  association  filter  (IMMJPDAF)  to  effectively  track  an  intercontinental  ballistic 
missile  (ICBM)  and  detect  decoys  and  maneuvers  is  the  focus  of  this  thesis.  Model 
development  and  data  association  schemes  are  discussed  along  with  optimized  values  for 
selected  parameters. 

Two  different  models  are  used  to  describe  the  missile  dynamics.  One  model 
describes  straight-line  motion  (SLM),  and  the  second  model  describes  a  constant  velocity 
turn.  The  interacting  multiple  model  (IMM)  algorithm  mixes  the  estimates  from  each  of 
the  models.  The  benefit  of  mixing  the  estimates  is  accurate  prediction  of  missile  position 
regardless  of  maneuvering  frequency.  Once  a  split  target  (e.g.,  decoy)  is  detected,  the 
common  history  between  the  missile  and  decoy  is  exploited  to  generate  a  separate  track 
for  the  decoy. 

The  joint  probabilistic  data  association  filter  (JPDAF)  matches  positions  from 
multiple  targets  received  by  the  sensor  to  the  predicted  location  of  each  existing  target. 
All  possible  combinations  of  received  positions  (observations)  and  predicted  target 
locations  are  assigned  a  value  based  on  probability  of  occurrence.  A  quantitative  analysis 
of  the  probability  values  allows  the  filter  to  distribute  observations  to  the  associated 
target.  The  successful  assignments  from  the  JPDAF  allow  the  IMM  to  better  predict  the 
next  target  position  and  increase  the  tracking  performance  of  the  integrated  IMMJPDAF. 

With  regard  to  ICBM  tracking,  the  IMM  algorithm  proved  to  be  capable  during 
both  steady-state  and  maneuvering  conditions.  For  the  given  scenarios,  the  estimates 
resulting  from  the  IMM  mixing  process  received  a  reduction  of  measurement  error  of  50 
percent  on  average.  The  degradation  of  the  individual  model  estimates  was  minimal.  The 
number  of  transitions  between  the  models,  especially  during  a  constant  velocity 


maneuver,  was  frequent  and  often  occurred  between  each  sample.  For  targets  with  the 
likelihood  of  transitioning  between  established  models  often,  this  method  is  appropriate. 
However,  for  the  purposes  of  an  ICBM  that  prefers  to  deploy  decoys  to  conducting 
maneuvers,  the  necessity  to  account  for  continuous  contact  maneuver  is  reduced.  In 
addition,  the  distance  errors  seen  during  model  transitions  are  concerning  due  to  the 
likelihood  of  decoy  deployment  during  those  periods.  Distance  errors  at  model  transitions 
had  negligible  effect  on  the  performance  of  the  IMMJPDA  filter  during  the  scenarios 
tested  but  could  become  a  concern  under  higher  split  target  densities. 

With  regard  to  split  target  tracking,  the  IMM  provided  a  convenient  framework  in 
which  to  generate  a  split  target  and  allowed  for  a  seamless  transition  for  follow-on  track 
maintenance.  We  introduced  a  non-traditional  IMM  approach,  where  the  two  models 
regarding  the  split  target  are  not  subject  the  standard  IMM  processes.  Initially,  stacking 
the  combined  estimates  to  generate  the  split  track  allowed  the  common  history  of  the 
target  and  decoy  to  be  exploited.  The  decoy  track  generation  caused  an  increase  in 
distance  error  of  the  combined  estimates  of  over  two  and  a  half  times  measurement  error. 
The  decoupling  of  model  three  and  model  four  with  the  target  models  prevented  further 
degradation  of  the  combined  estimates.  Additionally,  the  linearly  increasing  model  three 
errors  showed  that  earlier  decoupling  of  the  split  target  and  true  target  was  best. 

Performance  comparisons  of  the  resultant  algorithm  to  a  standard  Kalman  filter 
utilizing  a  nearest  neighbor  discriminator  are  conducted.  Scenarios  included 
combinations  of  missile  maneuver  and  decoy  deployment.  While  the  Kalman  filter 
experiences  limited  success,  the  IMMJPDAF  tracks  the  missile  in  every  scenario. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

Ballistic  missile  defense  is  not  a  new  concern  for  the  United  States,  but  until 
recently  the  number  of  countries  with  credible  threats  has  been  limited.  The  availability 
of  ballistic  missile  technology  to  unstable  nations  poses  a  legitimate  threat  to  the  security 
of  the  United  States.  The  challenge  of  hitting  high-speed  missiles  in  an  exo-atmospheric 
environment  has  been  overcome  with  capable  interceptors.  The  use  of  decoys  released 
from  the  missile  body  has  the  potential  to  dramatically  reduce  the  effectiveness  of  the 
current  interceptor  inventory.  Solutions  include  development  of  a  new  platform  or  an 
upgrade  to  existing  interceptor  sensors  and  tracking  algorithms.  The  more  fiscally 
prudent  solution  is  the  latter. 

In  response  to  growing  threats,  the  Ballistic  Missile  Defense  Agency  has  adopted 
a  multi-layered  defense  strategy  utilizing  the  full  spectrum  of  sensors  and  interceptors  for 
boost,  mid-course,  and  terminal  phase  interception  [1],  The  collection  of  sensors  and 
platforms  is  illustrated  in  Figure  1.  A  clear  distinction  between  the  phases  of  ballistic 
missile  travel  helps  to  divide  the  unique  challenges  posed  in  each.  The  boost  phase  is 
from  launch  to  the  end  of  powered  flight  and  lasts  between  60  and  300  seconds  [2], 
Interception  in  this  phase  is  preferred  but  challenging  based  on  the  short  window  of 
opportunity  and  political  issues  surrounding  interceptor  basing.  Mid-course  describes  the 
trajectory  from  booster  burnout  to  return  to  the  earth’s  atmosphere  [2],  Recognizing  the 
vulnerability  while  in  the  midcourse  phase,  threat  missiles  deploy  decoys  in  this  phase 
only  for  protection  during  the  100  second  to  10  minute  travel  time.  The  terminal  phase  is 
a  last  resort  for  interception  and  is  defined  where  gravitational  effects  are  noticed, 
typically  between  altitudes  of  70  and  100  km  [2], 
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Figure  1 .  A  combination  of  sensors  and  platforms  in  a  multi-layered  missile  defense 
strategy  contribute  to  successful  interception  of  a  ballistic  missile  in  each 

phase  of  travel.  From  [1], 


The  National  Research  Council  submitted  a  report  in  [2]  describing  the  current 
challenges  of  Ballistic  Missile  Defense.  Midcourse  intercept  was  identified  as  critical  to 
the  success  of  ballistic  missile  interception,  even  in  the  terminal  phase  [2],  If  destruction 
of  reentry  vehicles  that  split  off  from  the  missile  cannot  be  completed  in  the  mid-course 
phase,  identification  and  tracking  must  be  accomplished  for  terminal  phase  interception. 
Mid-course  discrimination  is  accomplished  by  the  Ground-Based  Interceptor  system 
utilizing  secure  satellite  communications  for  initial  threat  missile  parameter  evaluation 
and  information  dissemination  [1], 

B.  OBJECTIVE 

Advanced  missile  systems  are  capable  of  deploying  decoys  and  maneuvering  to 
thwart  midcourse  interception.  Algorithms  that  can  track  through  maneuvers  and  deal 

2 


with  a  single  target  splitting  into  multiple  independent  targets  when  decoys  are  deployed 
are  required  to  deal  with  these  advanced  threats.  The  feasibility  of  using  the  interacting 
multiple  model  (IMM)  [3]  and  the  joint  probabilistic  data  association  (JPDA)  [4]  tracking 
filters  in  a  ballistic  missile  interception  scenario  is  the  focus  of  this  research.  More 
specifically,  the  ability  to  accurately  predict  the  missile  state  during  maneuvering  and 
target  split  operations  will  be  analyzed. 

C.  PREVIOUS  WORK 

Ballistic  missile  tracking  algorithms  are  relied  upon  to  efficiently  process  the 
information  provided  by  the  vast  network  of  sensors  seen  in  Figure  1.  Effective 
interception  requires  accurate  threat  missile  position  and  speed  information.  Application 
of  the  probabilistic  data  association  filter  (PDAF)  in  over-the-horizon  radar  systems,  as 
discussed  in  [5],  provides  threat  missile  track  generation  from  a  cluttered  environment. 
To  increase  confidence  in  target  track  data,  the  integration  of  additional  target 
information  such  as  mass,  heat,  and  signal  strength  has  been  extensively  researched.  For 
example,  Colegrove  and  Davey  use  signal-to-noise  ratio  in  [6]  to  quickly  formulate  target 
tracks. 

The  IMM  best  fits  applications  involving  fast  moving  and  frequently  maneuvering 
targets  (e.g.,  air  traffic  control  systems,  GPS  navigation,  and  radar  tracking  systems). 
Ballistic  missile  tracking  using  the  IMM  has  been  accomplished  with  a  variety  of 
approaches  with  a  distinctive  feature  being  the  number  of  phases  included.  For  example, 
Cooperman  presents  a  model  for  all  three  phases  of  ballistic  travel  in  [7]  using  multiple 
sensors.  Farrell  uses  a  similar  approach  in  [8]  by  applying  standard  Kalman  filters  in  each 
of  his  separate  models  (constant  axial  force,  ballistic  acceleration,  and  autocorrected 
acceleration).  Both  approaches  achieved  position  error  results  that  are  best  suited  for  a 
single  phase  of  missile  travel,  while  the  other  phases  suffer  from  large  inaccuracies.  A 
majority  of  approaches,  such  as  the  one  found  in  [9],  use  multiple  models  to  describe  a 
single  phase,  and  most  often  the  focus  is  the  boost  or  terminal  phase.  Our  research  will 
focus  on  only  mid-course  missile  travel,  requiring  missile  models  with  reliable  split  target 
detection. 
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Methods  for  detection  of  split  targets  are  less  researched  due  to  limited 
applications.  A  method  for  split  target  tracking  in  clutter  using  an  IMM  is  presented  in 
[10],  but  the  assumptions  used  in  this  approach  limit  the  number  of  target  splits  (such  as 
weapons  shot  from  fighter  aircraft).  A  generalized  method  for  integrating  the  interacting 
multiple  model  joint  probabilistic  data  association  filter  (IMMJPDAF)  with  target  splits 
is  presented  in  [3]  and  is  the  basis  for  the  approach  used  in  this  research.  To  the  best  of 
the  author’s  knowledge,  the  application  and  parameter  optimization  of  the  IMMJPDAF 
for  tracking  a  maneuvering  ballistic  missile  during  decoy  deployment  has  not  been 
researched. 

D.  THESIS  ORGANIZATION 

The  standard  Kalman  filter,  IMM,  and  JPDA  equations  are  presented  and  adapted 
for  efficiency  in  Chapter  II,  along  with  methods  for  split  track  generation  and  deletion. 
The  generation  and  transition  between  models  are  discussed  in  Chapter  III  and  are  shaped 
by  an  estimation  of  the  most  likely  scenarios.  Critical  initiation  values  for  each  model  to 
include  sensor  and  plant  noise  characteristics  are  also  identified.  Performance  metrics  and 
an  analysis  of  MATLAB  simulation  results  are  given  in  Chapter  IV.  The  simulations 
presented  include  a  stand-alone  evaluation  of  the  IMM  during  target  maneuver,  split 
target  detection  from  various  target  states,  and  a  comparative  analysis  of  the  combined 
IMM  and  JPDAF  to  the  standard  Kalman  filter.  A  summary  of  results  and 
recommendations  for  future  work  are  the  focus  of  Chapter  V. 
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II.  FILTERING  METHODOLOGY 


A.  KALMAN  FILTER 

A  brief  overview  of  the  Kalman  filter  equations  is  given  here  and  taken  directly 
from  [11].  The  target  state  at  time  k  is  x(k),  representing  the  true  target  location.  The 
transition  of  the  true  state  of  the  system  is  described  by 

x(k  +  \)=F(k)x(k)+v(k),  (1) 

where  F(k )  is  the  state  transition  matrix  and  v(&)is  the  plant  noise.  The  deterministic 
term  u(k),  which  accounts  for  any  known  inputs  to  the  system  dynamics,  and  the 
associated  input  gain  matrix  G(k )  are  normally  included  in  Equation  (1)  but  are  not 
required  for  our  analysis.  Measurements  are  received  at  each  time  k  and  are  related  to 
the  true  target  state  as  follows: 

z(k+l)  =  H(k+l)x(k+l)+w(k+l),  (2) 

where  w(k  + 1)  represents  zero-mean,  white  measurement  noise.  The  state  estimate  at 
time  k+ 1  is 

x(k+l\k)  =  F(k)x(k\k),  (3) 

where  x(k  \  k )  is  the  previous  state  estimate.  Once  a  measurement  at  k  +  \  is  received, 
the  predicted  measurement  is  compared  to  the  received  measurement  as  follows: 

z(£+l)  =  z(£+l)  — z(£+l  |  k),  (4) 

where  the  predicted  measurement  is 

~z(k+l\k)  =  H(k+l)x(k  +  \\k).  (5) 

The  state  estimate  is  then  updated  by  incorporating  z(k  + 1)  in  the  following  manner: 

x(k+l\k+l)  =  x(k+l\k)  +  W(k+ \jz(k  + 1) ,  (6) 
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where  W(k+ 1)  is  the  Kalman  filter  gain.  The  Kalman  gain  requires  computations 
involving  the  state  covariance  P(k  \  k )  and  the  state  transition  matrices.  The  state 
prediction  covariance  is 

P(k+l\k)  =  F(k)P(k  |  k)F(k)'+  Q(k) ,  (7) 

where  Q(k)  is  the  covariance  of  the  zero-mean  white  process  noise.  Measurement  noise 
covariance  is  taken  into  consideration  as  follows: 

S(k  + 1)  =  H(k  + 1  )P{k + 1 1  k)H(k + 1)’+  R(k + 1) ,  (8) 

where  f?(£  +  l)  is  the  zero-mean,  white  measurement  noise  covariance.  The  Kalman  filter 
gain  calculation  is 

W(k  +  1)  =  .P(£+1|  k)H(k  + + 1) 1 .  (9) 

The  state  covariance  is  updated  using  W(k  + 1)  as  follows: 

P(k+l\k+l)  =  P(k+l\k)-W(k+l)S(k+l)W(k+l)'.  (10) 

Equation  (3)  through  Equation  (10)  are  performed  on  each  cycle. 

B.  INTERACTING  MULTIPLE  MODEL 
1.  Purpose 

Kalman  filter  tracking  errors  are  known  to  increase  as  a  target  continues  to 
deviate  from  the  trajectory  described  by  the  state  transition  matrix  shown  in  Equation  (1). 
Thus,  a  single  Kalman  filter  is  insufficient  for  tracking  a  maneuvering  target  or  multiple 
targets  regardless  of  their  trajectories.  The  IMM  approach  circumvents  these  problems  by 
creating  multiple  state  models.  A  model  is  defined  for  each  target  state  (e.g.,  straight-line 
motion,  constant  turn  rate)  to  best  describe  the  observed  motion  in  the  form  of  multiple 
transition  matrices.  Multiple  Kalman  filters  running  in  parallel  reduces  tracking  errors, 
assuming  an  adequate  algorithm  exists  for  filter  transitions.  “The  unique  feature  of  the 
IMM  approach  is  the  manner  in  which  the  state  estimates  and  the  covariance  matrices 
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from  these  multiple  models  are  combined  according  to  a  Markov  model  for  the  transition 
between  target  maneuver  states”  [4], 

The  IMM  flowchart  in  Figure  2  below  appears  in  [4]  and  provides  a  convenient 
method  of  summarizing  IMM  operation.  Once  a  measurement  (observation)  is  received, 
the  JPDAF  is  assigned  the  task  of  deciding  if  and  to  which  established  track  it  belongs, 
satisfying  the  “Gating  and  Data  Association”  block.  Explanation  of  IMM  operation 
begins  with  appropriately  assigned  observed  data  and  follows  the  flowchart  in  Figure  2. 


Model 

Probabilities 

- ► 


" Filtered  State  Estimates  and 
Covariance  Matrices 
X  j(k-1 1  k-1),  Pj(k-1|  k-1) 


IMM  Mixing 


Updated  Model  Probabilities, 
State  Estimates  and  Covariance 
Matrices 


State  and  Covariance 
Prediction 


l  X  j  (kl  k-1) ,  Pj(k  |  k-1 ) 


Gating  and  Data  Association 

Observation  Data 

^ - - 

|  Observation-to-T rack  Assignments 

Update  Model  Probabilities 

i 

Update  Filtered  Estimates 

- ^Output  to  User 

Figure  2.  An  IMM  flowchart  depicts  the  major  steps,  inputs,  and  outputs  of  the  IMM 

operation.  From  [4], 


2.  Model  Probabilities 

The  number  of  possible  target  states  determines  the  number  of  models  r  required 
for  the  IMM.  A  sufficient  collection  of  models  is  required  to  cover  the  spectrum  of 
anticipated  target  motion.  The  transition  between  states  is  described  by  both  conditional 
probabilities  and  transitional  probabilities.  Conditional  probabilities  ut  for  i  =  l,2,...,r 
are  initially  assigned  a  priori  values,  which  are  updated  at  each  step  Update  is  based 
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on  a  comparison  of  the  received  measurement  z{k)  to  each  of  the  respective  model  state 

estimates  Xt(k  |  k  —  1) .  The  probability  that  the  target  is  in  model  state  i  as  computed  just 
after  measurement  data  are  received  at  time  k—l  is  defined  as  ut{k— 1)  [4],  The 
normalized  statistical  distance  is  calculated  for  use  in  update  of  ut  and  is 

df=z(k)'S;\k)z(k),  (ii) 


where  the  residual  covariance  matrix  St(k)  =  HP(k\  k-  l)H'+R  is  received  from 
Equation  (8).  The  model  probability  update  utilizes  a  likelihood  function  based  on 
Gaussian  statistics  defined  as  follows: 


4(*)  = 


(12) 


where  M is  the  dimension  of  St(k) .  The  likelihood  function  is  used  to  update  the  model 
probabilities  by 


«(*)= 


4(*)C,(*) 
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(13) 


where  C  is  a  normalizing  constant  from  [4]  and  is  described  by 


1=1 


(14) 


The  variable  Ct(k- 1)  in  Equation  (14)  is  a  linear  combination  of  both  conditional  and 
transitional  probabilities  and  updated  during  the  IMM  mixing  process.  The  updated 
model  probabilities  from  Equation  (13)  assign  a  weighted  contribution  of  the  updated 
states  estimates  from  the  Kalman  filter  models.  Combined  estimates  are  generated  that 
represent  the  most  up  to  date  estimates  based  on  the  current  measurements.  The 
combined  estimates  and  associated  covariance  matrices  are  calculated  as  follows: 


xc(k\k-l)  =  YJu,(k-  1)*<  (k)  (15) 

1=1 


and 
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£<*)  =  £ii,W*|*-l)+[»(*|t-l)-x«(*-l)]Wi|*-l)-Jct(*-l)],)>  (16) 

i  1 

where  xc(k  —  1)  is  the  combined  estimate  from  the  previous  time  step,  x,(k  \  k  - 1)  is  the 
estimate  from  the  associated  model,  and  /*(&|  k—  1)  is  the  covariance  matrix  produced 
from  the  standard  Kalman  filter  equations  in  each  respective  model  [11],  At  this  point, 
the  model  probabilities  and  filtered  estimates  have  been  updated  based  on  observation 
data,  and  all  required  data  for  IMM  mixing  is  available. 

3.  IMM  Mixing 

The  IMM  mixing  satisfies  several  objectives.  First,  “the  IMM  approach  provides 
the  most  effective  framework  for  adaptive  filtering.”[4]  Each  set  of  observation  data 
allows  discrimination  or  favoring  of  the  model  that  it  is  most  closely  associated  with. 
Additionally,  “the  mixing  process  uses  the  accurate  filter  state  estimates  to  correct  the 
less  accurate  filter  state  estimates”  [4],  The  process  analyzes  each  possible  target 
transition  and  adjusts  the  state  estimates  based  on  the  most  likely  target  state.  Mixing 
occurs  after  update  of  the  state  estimate  and  covariance  matrices  but  before  prediction. 

A  Markov  transition  matrix  is  utilized  to  describe  the  probability  pi}  that  “the 
target  will  make  the  transition  from  model  state  /  to  model  state  j  ”  [4],  For  a  target 
characterized  by  three  states,  the  transition  matrix  is 


Pu 

Pl2 

Pl3 

P  21 

P  22 

P'23 

Pi  1 

P  32 

P  33 

The  probabilities  contained  within  M  are  developed  a  priori  and  remain  constant.  Each 
row  must  sum  to  unity  to  prevent  degradation  of  the  estimates  contributed  by  the 
associated  model. 

While  the  Markov  matrix  describes  the  estimated  transitions  between  models,  a 
parameter  describing  actual  recent  model  transitions  is  required  to  update  model 
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probabilities.  “The  conditional  probability  that  the  target  made  the  transition  from  state  / 
to  state  j  at  time  A: —  1  is  utJ{k  —  1)’  [4].  Conditional  probability  differs  from  py  in  that 

the  transition  from  state  /  to  state  J  occurs  given  that  the  target  is  initially  in  state  i. 
The  relationship  between  the  conditional  and  transitional  probability  is 


(*-!)  = 


C,{k- 1)  ’ 


(18) 


where  Cj(k— 1)  is  the  probability  after  interaction  that  the  target  is  in  state  j  and  can  be 
defined  as 


C,<i-l)  =  2>#A(i-l).  (19) 

i=l 

Once  transitional  probabilities  are  calculated  for  each  state,  a  mixing  process  is 
executed  to  produce  new  filtered  state  estimates  and  Kalman  filter  covariance  matrices. 
For  all  models  j  =  1,2,. ..,r  the  state  estimates  and  covariance  matrices  are  calculated 
using  intermediate  values  described  by 

n 

Xj  {k  - 1 1  k  - 1)  =  ^ py(k  -  1)jc,-(£  - 1 1  k  - 1)  (20) 

il 

and 

1=1 

The  IMM  mixing  is  complete. 

The  prediction  step  used  in  the  Kalman  filter  is  then  applied  to  produce  a 
predicted  state  estimate  and  a  Kalman  filter  covariance  matrix  as  follows: 

xi(k+l\k)  =  Fixt  (22) 

and 

Pi(k  +  l\k)  =  FiPI°Fi'+Qr  (23) 
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The  “State  and  Covariance  Prediction”  block  of  the  IMM  flowchart  in  Figure  2  is 
satisfied  by  Equation  (22)  and  Equation  (23).  Upon  receipt  of  new  observation  data,  the 
IMM  process  is  repeated. 

C.  JPDAF 

Ideally,  no  ambiguity  exists  in  regards  to  which  track  is  associated  with  a  given 
observation;  however,  such  an  assumption  is  an  oversimplification  in  tracking  scenarios. 
Proper  assignment  of  measurement  data  to  the  correct  target  is  vital  to  maximizing  the 
performance  of  both  the  Kalman  filter  and  the  IMM.  To  address  issues  of  track 
assignment,  the  probabilistic  data  association  (PDA)  method  associates  measurement  data 
with  established  tracks  and  removes  extraneous  measurements.  The  JPDAF  is  an 
extension  of  the  PDA  to  a  multiple  target  scenario  and  serves  to  resolve  any  existing 
correlation  conflict.  JPDAF  integration  with  the  IMM  is  described  in  Section  D. 

We  will  use  an  example  posed  in  [4]  and  illustrated  in  Figure  3  below  to  describe 
the  operation  of  the  JPDAF.  The  scenario  assumes  two  tracks  have  been  established.  The 
variables  PI  and  P2  in  Figure  3  represent  predicted  states  from  the  two  established  tracks. 
Both  predictions  are  located  at  the  center  of  the  respective  circular  gates  and  are 
established  along  each  target’s  associated  track.  A  circular  gate  is  used  in  this  example 
for  ease  of  computation.  Gate  size  is  dependent  on  noise  characteristics  and  potential 
target  maneuver  parameters  defined  within  the  target  model.  Three  observations  are 
received  simultaneously  and  are  denoted  by  01,  02,  and  03.  All  three  observations  are 
considered  for  validation  since  each  is  located  within  one  or  both  track  gate(s).  The 
probability  that  a  true  target  observation  will  fall  within  the  correct  track  gate  PG  is 
assumed  to  be  one. 
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Track  gates 


Figure  3.  A  typical  JPDAF  scenario  shows  the  need  for  conflict  resolution.  PI  and  P2 
are  the  predicted  track  positions.  01,  02,  and  03  are  the  three  observations 
that  fall  within  the  track  gates.  From  [4], 

A  hypothesis  matrix  is  constructed  to  track  each  possible  observation  and  track 
combination,  including  the  possibility  that  none  of  the  observations  belong  to  either  track. 
A  probability  is  calculated  for  each  hypothesis  based  on  a  Gaussian  likelihood  function. 
Each  track  i  is  related  to  each  observation  J  by 


d.n 
e  5 


(24) 


where  M  is  the  measurement  dimension  [11],  To  fully  describe  the  probability  of  the 
hypothesis,  the  probability  of  detection  PD  and  a  false  return  density  function  fi  are 
included.  The  full  development  of  the  probability  equations  are  given  in  [4],  but  the 
equations  used  for  the  scenarios  dealt  with  in  this  research  are  contained  in  Table  1. 
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Table  1.  A  JPDAF  hypothesis  matrix  is  constructed  for  the  scenario  presented  in 
Figure  3  with  a  probability  of  detection  equal  to  one.  From  [4], 


Hypothesis  Number 

Track  Number 

Hypothesis  Likelihood 

PI 

P2 

h\Ht) 

1 

01 

02 

S\\S22^dP 

2 

03 

02 

S\?,S22^dP 

3 

01 

03 

&\\&2?> PdP 

4 

02 

03 

&\2&2?PdP 

The  hypotheses  in  Table  1  are  constructed  by  listing  all  possible  outcomes  under 
each  track.  The  scenario  in  Figure  3  shows  that  all  three  observations  (01,  02,  and  03) 
may  be  associated  with  track  one  (PI),  and  only  observations  two  and  three  may  be 
associated  with  track  two  (P2).  Observation  one  does  not  appear  in  the  track  two  column, 
but  all  other  combinations  of  tracks  and  observations  are  represented.  For  example, 
hypothesis  one  states  that  observation  one  (01)  is  assigned  to  PI  and  that  observation  two 
(02)  is  assigned  to  track  two  (P2).  Observation  one  lies  outside  the  gate  of  P2  and  is, 
therefore,  not  assigned  to  P2.  We  will  assume  that  the  probability  of  detection  PD  is  one, 
meaning  that  each  track  always  receives  a  measurement.  This  is  a  significant  assumption 
and  greatly  reduces  the  number  of  hypotheses  that  must  be  considered.  For  space-based 
sensor  tracking  of  an  exo-atmospheric  ballistic  missile,  the  assumed  P  is  not  unrealistic. 


The  hypotheses  likelihoods  are  computed  but  require  normalizing.  Normalized 
probabilities  are  calculated  by 


m  ? 

2 >w,) 


J=1 


(25) 


where  m  is  the  total  number  of  hypotheses.  The  probability  that  01  is  associated  with  PI 
is  calculated  using  Table  1  as  follows: 
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(26) 


where  /?,,  is  the  probability  that  01  belongs  to  PI.  The  remaining  probabilities  for  track 
one  ( hu  and  /zB)  are  computed  in  the  same  manner  and  must  all  sum  to  one.  Track  two 
probabilities  are  computed  the  same  way.  Notice  that  hn  is  zero  because  01  does  not  fall 
within  the  gate  of  P2. 

The  information  obtained  from  each  track’s  association  with  each  observation  is 
incorporated  into  the  standard  Kalman  filter  equations.  To  account  for  multiple 
observations  in  Equation  (4),  a  weighted  sum  of  residuals  is  taken  by 


N 


(27) 


where  zy(k)  =  z^k)-  Hxt{k  \  k  —  \).  The  state  estimate  in  Equation  (5)  is  calculated  by 

setting  Zi(k)  equal  to  v(£  +  1)  .  The  Kalman  covariance  matrix  requires  inclusion  of 
uncertainties  from  each  of  the  observations.  If  we  assume  a  valid  observation  falls  within 
each  track  gate  during  each  discrete  time  step,  the  Kalman  covariance  is  the  sum  of 
Equation  (9)  with 


N 


(28) 


D.  JPDAF  AND  IMM  COMBINATION 

The  JPDAF  and  the  IMM  integration  require  time  phasing  of  their  respective 
parameters.  Model  probabilities  fiy(k— 1),  state  predictions  Xi(k  \  k  —  1) ,  and  Kalman 
filter  covariance  matrices  F*(/t|£— 1)  are  required  from  the  IMM  during  each  discrete 
time  step.  The  JPDAF  is  responsible  for  determining  the  gate  regions  for  each  model  and 
validating  observations  for  track  update  that  fall  in  each  model  gate.  Blackman  and 
Papoli  outline  a  process  discussed  below  for  integration  of  all  parameters  in  [4],  Once 
new  observations  are  received,  the  IMM  submits  a  likelihood  function  for  each  model: 
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(29) 


N  p  -4n 

4=d  -Pd)J3+X-  d 


where  N  is  the  total  number  of  observations.  The  inclusion  of  all  hypotheses  defined  in 
the  JPDAF  causes  Equation  (24)  to  differ  from  the  previously  defined  likelihood  function 
in  Equation  (16).  Model  probabilities  are  updated  using  Equation  (17).  The  state 
estimates  and  the  covariance  matrices  are  then  updated  for  each  model  as  discussed  in  the 
JPDAF  following  Equation  (22).  The  final  step  includes  the  final  steps  of  the  IMM 
mixing  process  to  produce  new  filtered  state  estimates  and  Kalman  filter  covariance 
matrices  as  described  in  Equation  (14)  and  Equation  (15). 


E.  SPLIT  TARGET  TRACK  GENERATION 

The  example  used  above  to  describe  the  JPDAF  operation  dealt  with  track 
maintenance  and  ignored  the  formation  of  tracks.  We  will  assume  that  the  missile’s  track 
is  reasonably  developed  based  on  external  radar  information  and  that  other  targets  of 
interest  are  in  the  vicinity  of  the  missile  (e.g.,  deployed  decoys).  Thus,  any  measurements 
that  occur  within  the  missile’s  gate  region  and  that  cannot  be  associated  with  an  existing 
track  will  be  a  candidate  for  a  split  target.  An  automatic  track  formation  process  is 
described  in  [10],  but  we  will  apply  a  split  model  within  the  IMM  framework  to  support 
track  generation. 

The  modeling  of  decoy  deployment  follows  the  methodology  described  in  [3].  For 
our  purposes,  a  decoy  is  synonymous  with  a  split  target.  Tracking  a  split  target  adds  two 
models  to  the  IMM.  The  first  model,  the  “just  split”  model,  contains  a  stacked  state 
vector  and  takes  advantage  of  the  cross-covariance  terms  containing  the  past  common 
history  of  the  two  targets.  The  second  model  associated  with  the  split  continues  track 
maintenance  in  the  same  manner  that  the  missile  track  is  being  updated.  All  splits  are 
contained  within  this  model,  and  the  JPDA  is  used  to  discriminate  between  split  targets. 

The  “just  split”  model  described  below  is  discussed  in  [3]  and  contains  several 
assumptions.  The  model  begins  with  a  single  predicted  measurement  from  the  previous 
time  step,  and  two  validated  measurements,  as  discussed  in  the  JPDA,  are  present  within 

15 


the  missile  target  gate  at  the  current  time  step.  Scenarios  containing  greater  than  two 
measurements  (even  if  two  splits  are  occurring  simultaneously)  are  not  covered  in  this 
research.  Additionally,  we  will  assume  that  once  a  target  splits,  no  recombination  occurs; 
a  reasonable  assumption  due  to  the  likely  maneuver  of  the  missile  once  a  decoy  is 
deployed.  Therefore,  the  recombination  model  discussed  in  [3]  has  been  ignored. 

Once  two  measurements  are  received,  a  stacked  vector  consisting  of  a  duplicated 
state  is  created,  resulting  in 


xJ(k) = 


XJx{k) 

*2  (*). 


(30) 


where  J  represents  the  “just  split”  model  [3].  In  addition  to  the  same  state,  the  predicted 
state  is  identical  as  well  and  is  stacked  in  the  same  manner.  The  common  predicted 
measurement  then  becomes 


x\k  |  k)  =  x{k  |  jfc  - 1) + WJ(k)vJ(k) , 


(31) 


-J 


-J 


where  vJ(k)  =  zt(k)-z  (k),  z  (k)  =  HJ(k)x  (k\k-l),  and  the  just  split  Kalman  gain  is 


WJ(k)  = 


W/(k)  0 

0  W’{k) 


(32) 


Each  of  the  parameters  contained  within  the  Kalman  filter  equations  require  a 
similar  stacking  for  the  dimensions  of  the  matrix  calculations,  including  the  measurement 
noise  covariance  RJ(k )  and  the  measurement  matrix  FJ(k).  The  covariance  associated 
with  the  prediction  will  be  constructed  by 

P,i(k \k-\)  F^(k\k- 1)  1 

P'(*|*-l)  Pi(k\k-l)  J' 

The  combined  covariance  matrix  from  Equation  (16)  serves  as  the  initial  value  for  each 
of  the  covariance  terms  listed  in  Equation  (33).  The  off-diagonal  covariance  matrices 

contain  the  correlation  between  the  new  estimates  and  the  past  common  history  [3], 
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PJ(k\k-l)  = 


With  the  formation  of  the  split  track  underway,  we  are  now  interested  in 
decoupling  the  targets  and  tracking  each  independently.  The  validation  regions  discussed 
in  the  JPDA  become  useful  metrics  in  determining  a  “no  overlap”  test  explained  in  [3]. 
The  already  developed  missile  gate  region  is  compared  to  a  newly  developed  split  target 
gate  region.  As  the  targets  separate  such  that  the  two  gate  regions  no  longer  overlap,  we 
assume  the  two  targets  are  completely  decoupled.  The  “no  overlap  test  can  be  described 
by 

[^f  (*  |  *  - 1)  -  |  *  - 1)]'[^(^)  +  ^2(*)]  '[zf  (*  |  *  - 1)  -  (*  |  *  - 1)]  >  (34) 

where  g  is  chosen  as  the  tail  end  of  a  chi-squared  density  such  as  one  percent  [3].  Once 
the  threshold  for  decoupling  is  met,  both  targets  can  be  tracked  independently  with 
existing  state  estimates  and  covariance  matrices. 

F.  TARGET  DELETION 

We  have  chosen  to  suspend  split  target  tracking  at  a  specified  distance  between 
the  split  target  and  the  missile.  The  split  target  is  likely  to  transit  outside  the  interceptor 
field-of-view  before  software  interaction  is  required.  However,  split  target  tracks  may  be 
of  interest  for  kill  vehicle  adjudication  or  to  simplify  the  JPDA  problem  should  the  split 
target  return  within  the  missile  gate.  For  the  purpose  of  the  simulation,  we  choose  a 
threshold  range  from  the  missile  in  which  to  discontinue  tracking  of  80,000  meters; 
although,  split  target  measurements  continue  to  be  plotted. 
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III.  MODEL  CONSTRUCTION 


A.  OVERVIEW 

One  master  script  file  and  several  function  files  were  constructed  and  run  in 
MATLAB  2012a.  We  assume  measurement  data  is  converted  from  likely  range  and 
bearing  data  to  a  Cartesian  coordinate  system;  we  ignore  mathematical  errors  associated 
with  such  a  conversion.  A  variable  naming  scheme  is  used  throughout  the  MATLAB  files 
because  the  same  parameters  are  duplicated  between  models.  In  general,  a  number 
corresponding  to  an  associated  model  follows  each  repeated  variable  (i.e.,  h\  denotes  the 
state  transition  matrix  for  model  three).  All  units  are  in  meters. 

We  chose  four  models  to  fully  describe  missile  and  split  target  tracking. 
Assumptions  include  point  mass  modeling  of  both  the  target  missile  and  the  split  target 
with  motion  described  by  linear  control  dynamics.  Model  one  consists  of  a  missile  in 
straight-line  motion  at  constant  velocity.  Model  two  simulates  missile  maneuver  at  a 
constant  turn  rate.  Model  three  detects  and  analyzes  for  a  potential  target  split  from  the 
missile.  Lastly,  model  four  maintains  track  on  all  confirmed  target  splits  and  discontinues 
tracking  when  the  80,000-meter  threshold  is  reached. 

B.  MASTER  SCRIPT  FILE 

The  master  script  file  defines  common  variables  utilized  by  each  of  the  models, 
generates  true  target  vectors,  and  calls  function  files  as  required.  The  initialization  values 
chosen  are  discussed  as  we  present  them  and  are  embedded  throughout  the  model 
construction  explanations.  We  used  a  sample  time  of  0.01  seconds  over  1,000  samples, 
resulting  in  a  total  simulation  time  of  10  seconds. 

1.  True  Missile  Motion 

The  missile  vector  is  modeled  as  a  point  mass  containing  two  degrees  of  freedom.  At  the 
beginning  of  each  simulation,  the  true  target  state  vector  is  set  to 
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X 

0 

Vx 

2500 
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0 

Vy 

2300 

z 

1000000 

Vz 

0 

(35) 


where  x,y,z  describe  the  position  in  the  Cartesian  plane  and  vx,v  ,vz  describe  the 

velocity  in  each  dimension.  The  target  vector  is  initialized  to  resemble  a  DF-41  Chinese 
ICBM  with  a  velocity  of  2670  m/s  and  a  constant  altitude  of  1,000  km  and  is  taken  from 
open  source  information.  The  state  transition  matrix  for  straight-line  motion  is 


1  delta  0  0  0  0 

0  1  0  0  0  0 

0  0  1  delta  0  0 

0  0  0  1  0  0 

0  0  0  0  1  delta 

0  0  0  0  0  1 


(36) 


where  delta  is  the  discrete  step  time  of  0.01  seconds  [12].  Thus,  SLM  results  as  each 
time  step  increments  the  x  position  vector  by  25  meters  and  the  y  position  vector  by  23 
meters.  Each  six  by  one  true  target  motion  vectors  for  SLM  is  collected  in  a  matrix 
throughout  the  simulation.  Each  column  represents  the  target  state  at  a  specified 
measurement  time  step. 

Although  the  simulation  models  motion  and  measurements  in  all  three  dimensions, 
the  scenarios  analyzed  below  assume  a  constant  value  of  z  for  true  target  motion.  All 
plots  are  presented  in  the  x—y  coordinate  plane.  The  model  for  target  acceleration 
assumes  a  constant  speed  turn  in  the  x—y  plane  with  a  state  transition  matrix  as  follows: 
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(37) 


where  a  is  a  constant  turn  rate  and  b  is  the  discrete  time  step  [12],  The  turn  rate  value 
was  selected  such  that  a  20  g  turn  is  simulated.  The  turn  is  initiated  and  terminated 
according  to  specified  discrete  time  step  values.  Any  combination  of  SLM  and  turning  in 
either  direction  is  possible. 


2.  True  Split  Target  Motion 

Actual  decoy  deployment  is  triggered  at  a  specified  time  step  n  as  a  target  split. 
The  decoy  (split  target)  is  initialized  using  the  current  value  of  the  true  target.  The 
velocity  component(s)  of  the  decoy  are  then  modified  to  create  separation  from  the  true 
target  over  time.  We  assume  linear  motion  of  the  split  targets,  leading  to  the  same 
transition  matrix  for  F3  as  was  used  for  SLM  ( Fx ).  No  computational  changes  are  made 
as  the  true  split  vector  is  analyzed  in  model  four,  but  the  same  transition  matrix  is 
relabeled  as  F4  and  the  true  split  state  is  relabeled  as  x4  for  state  estimate  comparison 
purposes. 


3.  Sensor  Noise 

We  assumed  a  fixed  sensor  with  a  five-meter  error  in  each  dimension  (x,  y,  and 
z).  A  white  Gaussian  process  is  used  to  generate  the  sensor  noise  w(k)  from  Equation 
(2).  The  resultant  measurement  covariance  matrix  is 
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(38) 


~  2  2 

where  <r2  ,  cry  ,  and  <xz  are  equal  to  25  meters.  Each  of  the  Kalman  filters  operating  in 

parallel  utilizes  this  parameter.  “Typical  tracking  error  decreases  are  moderate  (on  the 
order  of  25  percent)... unless  very  accurate,  high  update  rate  measurement  data  are 
available.”[9]  Based  on  greater  than  a  40  percent  drop  in  mean  distance  errors  with  the 
specified  tracking  error  and  sampling  rate,  we  assess  our  estimations  to  be  reasonable. 


4.  Plant  Noise 

Plant  noise  v(k )  from  Equation  (1)  is  modeled  as  a  white  Gaussian  process  that  is 
independent  of  sensor  noise  and  is  used  to  describe  un-modeled  parameters  [12].  The 
plant  noise  covariance  matrix  is 
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(39) 


where  delta  is  the  discrete  time  step  of  0.01  seconds  per  sample  and  q2  is  a  scalar  value 
used  to  differentiate  between  the  different  model  states  [12].  The  value  of  q 2  is 
maintained  at  10  for  SLM  and  4,000  for  a  constant  speed  maneuver. 
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5. 


Model  Transitions 


We  restrict  the  transition  between  certain  models  based  on  both  physical 
constraints  and  assumptions.  We  first  assume  that  from  SLM  (model  one)  the  missile  is 
able  to  remain  in  SLM  with  probability pu,  conduct  a  maneuver  with  probability pvl,  or 
split  with  probability  pn  .  A  distinction  must  be  made  concerning  the  “transition”  from 
model  one  to  model  three.  We  are  not  saying  that  the  target  either  proceeds  in  SLM  or 
splits,  but  more  accurately,  that  a  split  may  occur  while  the  missile  proceeds  down  its 
track  without  maneuvering.  The  missile  is  physically  capable  of  doing  both,  and  missile 
state  estimates  continue  to  be  calculated  despite  a  split  occurring. 

We  will  also  assume  that  from  model  two  the  missile  is  able  to  return  to  SLM 
with  probability  p2l,  continue  to  maneuver  with  probability  p22,  or  conduct  a  target  split 
with  probability  pa .  The  same  argument  as  above  applies,  and  a  potential  target  split 
(model  three)  occurs  simultaneously  with  continued  missile  maneuver  (model  two).  We 
considered  setting  pn  to  zero  considering  the  most  likely  evasion  technique  to  be  a 
deployed  decoy  followed  by  a  missile  maneuver.  However,  an  equally  feasible  tactic  may 
include  deployment  of  a  decoy  during  a  constant  turn.  The  approach  used  by  Bar-Shalom, 
Chang,  and  Blom  in  [3]  only  account  for  target  splits  from  SLM  and  not  for  a  turn. 

Assuming  no  recombination,  transition  probabilities  for  model  three  and  four  are 
straightforward.  We  assume  potentially  split  targets  in  model  three,  allowing  only  options 
for  the  target  to  be  a  confirmed  split  ( pM )  or  just  an  erroneous  measurement  ( p33).  The 
transitional  probability  p33  also  accounts  for  track  formation,  as  a  finite  number  of 
measurements  are  required  to  confirm  a  split  target.  A  split  track  will  remain  in  model 
four  ( pu  )  for  track  maintenance,  resulting  in  a  transitional  probability  of  one. 

Now  that  we  have  identified  the  non-zero  transitional  probabilities,  values  must 
be  assigned  to  each.  A  lack  of  historical  ICBM  flight  and  decoy  data,  at  least  in  open 
sources,  forces  estimated  values.  Assuming  a  dominant  straight-line  trajectory  and  low 
split  target  density,  the  Markov  chain  transition  matrix  was  established  as 
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(40) 


0.8  0.1  0.1  0 

0.7  0.2  0.1  0 

0  0  0.3  0.7 

0  0  0  1 


An  unnecessary  amount  of  complexity  and  restriction  was  realized  soon  after 
model  construction  when  following  the  approach  presented  in  [3],  which  led  to  the 
generation  of  Equation  (40).  The  apparent  sacrifice  of  true  target  accuracy  as  the  number 
of  splits  increased  would  pose  an  issue  for  future  work.  Although  only  designed  to  handle 
a  limited  number  of  target  splits  in  this  research,  a  realistic  scenario  involves  frequent 
splits  or  batch  splits  over  the  duration  of  midcourse  flight.  A  rising  split  probability  ( pn 
and  )  results  in  a  reduction  of  other  model  transition  probabilities.  Instead  we  decided 
to  limit  the  IMM  model  transitions  to  model  one  and  model  two.  Multiple  observations 
falling  within  the  true  target  gate  will  now  trigger  a  transition  to  model  three,  a  potential 
target  split.  Thus,  the  new  transition  matrix  for  the  IMM  becomes 


Py  = 


0.8  0.2 
0.6  0.4 


(41) 


The  “modelprob”  function  file  utilizes  the  values  in  Equation  (40)  to  accomplish  IMM 
mixing. 


C.  MODEL  FUNCTIONS  AND  INTERACTIONS 

Each  model  generates  estimates  for  its  associated  trajectory.  The  process  used  to 
generate  estimates  from  the  true  target  vectors  is  described  below  for  each  track.  The 
master  script  file  triggers  when  each  model  (contained  within  function  files)  is  required 
based  on  the  discrete  time  step  k. 


1.  Straight-Line  Motion 

A  standard  Kalman  filter  is  applied  to  track  true  target  SLM.  Sensor  noise  is 
added  to  the  true  position  vector  by  multiplying  each  dimension  of  the  sensor  covariance 
by  a  pseudorandom,  normally  distributed  vector  of  the  same  size,  creating  measurements 
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with  simulated  zero  mean,  white  noise.  The  state  and  covariance  are  updated  in 
accordance  with  Equation  (6)  and  Equation  (10),  respectively.  A  distance  error 
calculation  is  performed  following  state  estimate  update.  The  model  outputs  the  state 
estimate  and  covariance,  the  noisy  measurements  for  use  in  model  two,  measurement 
error,  and  the  distance  error.  The  function  file  “modelprob”  collects  all  necessary  data 
from  model  one  and  model  two  to  conduct  IMM  mixing,  update  model  probabilities,  and 
update  the  fdtered  estimates.  The  intermediate  values  for  state  estimate  and  covariance 
presented  in  Equation  (20)  and  Equation  (21)  are  sent  out  from  “modelprob”  to  a  separate 
function  file,  “prediction”,  which  generates  the  final  model  estimates  from  Equation  (22) 
and  Equation  (23). 

2.  Constant  Speed 

The  second  model  receives  the  noisy  measurements  generated  in  model  one  and 
conducts  the  Kalman  filter  update  equations.  Distance  error  is  also  calculated  and  sent  as 
an  output,  along  with  the  updated  state  and  covariance  estimates. 

3.  Potential  Split  Target 

Model  three  requires  both  true  target  and  split  target  vectors  to  generate  noisy 
measurements  for  both.  The  model  one  function  file  is  not  called  when  analyzing  a  split. 
The  noisy  measurements,  along  with  other  required  variables,  are  stacked  and  prepared 
for  split  target  evaluation  as  discussed  in  Equation  (30)  through  Equation  (33).  A  JPDA 
calculation  is  conducted  utilizing  Equation  (24)  for  a  two  observations  and  two  track 
scenario.  The  measurement  vector  is  stacked  according  to  the  JPDA  filter  results,  and 
state  and  covariance  updates  are  calculated  while  stacked.  The  state  estimates  and 
associated  covariance  matrices  are  separated  out.  State  estimates  and  covariance  matrix 
predictions  are  calculated  for  model  three  and  are  available  for  output.  The  updated 
model  one  estimates  and  covariance  matrix  are  sent  to  “modelprob”  and  “prediction”  just 
as  when  no  split  is  being  analyzed  for. 
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4.  Split  Target  Tracking 

The  fourth  and  final  model  receives  the  un-stacked  state  estimates  from  model 
three  and  continues  track  maintenance  on  all  split  targets.  The  true  split  target  tracks  are 
received,  and  sensor  noise  is  added  to  the  measurements.  Standard  Kalman  filtering 
occurs  and  predicted  state  estimates  result.  When  the  received  measurement  exceeds  the 
state  estimate  by  greater  than  80,000  meters,  tracking  is  terminated. 
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IV.  MATLAB  SIMULATIONS 


A.  SIMULATION  SETUP 

Simulations  were  designed  to  test  the  ability  of  the  multiple  filter  models  to  track 
a  target  through  maneuvering  turns  and  maintain  the  correct  track  as  decoys  are  deployed. 
Mean  distance  errors  are  calculated  over  100  Monte  Carlo  simulation  runs. 

B.  INDEPENDENT  EVALUATIONS 
1.  IMM  Missile  Tracking 

Accurate  tracking  of  the  ICBM  through  a  turn  is  the  goal  of  the  first  simulation. 
The  scenario  establishes  a  true  target  track  beginning  in  SLM,  executing  a  turn  from 
sample  k  —  200  until  sample  k  —  800  and  then  returning  to  SLM.  The  combined  estimate 
track  is  used  as  a  measure  of  IMM  tracking  performance.  To  validate  the  additional 
computations  in  the  IMM  mixing  process,  the  individual  estimates  from  model  one  and 
model  two  are  compared  to  the  combined  estimates.  The  model  one  and  model  two 

estimates  x(&  +  l|£+l)  are  collected  following  the  Kalman  filter  update  (Equation  5), 

but  the  state  estimates  x(k + 1 1  k)  being  updated  were  filtered  from  the  IMM  mixing 
process  in  the  previous  time  step  k-  The  model  one  and  model  two  estimates  represent 
the  quantities  at  the  top  of  the  IMM  flowchart  in  Figure  2,  and  the  combined  estimates 
represent  the  updated  filtered  estimates. 

For  the  first  simulation,  the  number  of  samples  is  10,000  with  a  sample  rate  of 
100  samples  per  second.  The  true  target  motion  is  shown  in  Figure  4  on  an  r-y 
coordinate  plane  and  provided  as  an  overview  of  the  simulation.  The  boxes  in  Figure  4 
correspond  to  the  location  in  which  selected  focus  areas  in  the  scenario  are  taken  from. 
The  locations  include  steady  state  SLM,  the  onset  of  target  maneuver,  and  steady-state 
maneuver.  The  track  layouts  and  distance  error  plots  are  presented  and  discussed  at  each 
focus  area. 
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Figure  4.  True  target  motion  on  an  x-y  plane  identifying  the  location  of  areas  of  focus 

for  follow-on  figures. 


The  combined  estimates  provided  better  tracking  performance  in  both  straight- 
line  and  accelerating  scenarios  than  the  individual  state  estimates.  The  track  layouts 
displayed  in  Figure  5  contain  only  a  snapshot  of  the  relationship  between  the  tracks 
during  SLM.  The  position  of  each  track  relative  to  the  true  target  (the  ICBM)  track  is  as 
expected.  During  straight-line  motion,  model  one  is  expected  to  be  closer  to  the  true 
target  and  model  two  further  away.  The  combined  estimate  track  sits  closer  than  even  the 
straight-line  estimate  track  for  a  majority  of  the  time.  All  tracks  demonstrate  an 
improvement  over  the  received  measurements. 
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Figure  5.  True  target  track,  IMM  estimate  tracks,  and  measurements  during  straight-line 

motion  on  a  Cartesian  plane. 


The  track  layout  displayed  in  Figure  6  follows  a  similar  discussion  to  those  seen 
in  Figure  5  except  for  the  maneuvering  scenario.  Although  simply  a  snapshot  (the  layout 
between  two  measurements  is  shown),  the  tracks  maintain  this  pattern  for  a  majority  of 
the  time.  In  this  case,  model  two  estimates  are  closer  to  the  true  target  track,  followed  by 
the  combined  estimates,  then  track  one  estimates,  and  finally  the  measurements.  The 
combined  estimates  also  tend  to  lie  closer  to  the  model  one  estimates  despite  a  target 
maneuver  based  on  the  probability  assignments  given  in  Equation  (40). 
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Figure  6.  True  target  track,  IMM  estimate  tracks,  and  measurements  during  a  constant 

speed  maneuver  in  the  Cartesian  plane. 


The  results  contained  in  Figure  7  are  provided  to  capture  the  relationship  between 
the  estimates  and  further  prove  the  effectiveness  of  the  IMM  process.  The  tracks  progress 
from  the  upper  left  comer  of  the  figure  to  the  lower  right  portion.  As  expected,  the  model 
one  track  continues  to  be  furthest  from  the  true  target  track  during  a  constant  speed  turn. 
In  the  upper  left  portion  of  Figure  7,  the  measurements  are  near  the  true  track  and 
increase  the  model  two  probability  u2  from  Equation  (17).  The  values  contained  in 
Equation  (19)  and  in  Equation  (20)  are  driven  up  as  a  result  of  u2  going  up  and  cause  an 
immediate  shift  of  the  combined  estimates  towards  the  model  two  track.  As  the 
measurement  vector  diverges  from  the  true  track  and  toward  the  model  one  estimates,  u2 
goes  down  and  the  model  one  probability  ux  goes  up.  The  contribution  from  model  one, 
which  is  further  from  the  true  track,  drives  the  combined  estimates  away  from  the  true 
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track.  In  this  snapshot,  the  measurements  drive  the  location  of  the  combined  estimates 
within  a  channel  consisting  of  model  one  and  model  two.  This  process  is  repeated  on 
each  time  step  and  causes  slight  variations  in  the  combined  estimates. 


Figure  7.  True  target  track,  IMM  estimate  tracks,  and  measurements  during  constant 

speed  turn  in  a  Cartesian  plane. 


The  relationship  between  the  track  estimates  and  true  target  track  are  quantified 
by  distance  error.  The  summation  and  averaging  of  the  combined  estimate  variations  over 
several  (100  Monte  Carlo)  simulations  reveals  their  performance.  During  SLM  the  results 
in  Figure  8  are  as  expected.  Model  one  contains  the  least  error,  followed  by  the  combined 
estimates,  and  then  model  two.  Estimated  distance  error  values  exhibiting  a  reduction  of 
nearly  50  percent  of  measurement  error  validate  performance  of  the  IMM  filter. 
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Figure  8.  SLM  mean  distance  errors  of  IMM  estimates  with  mean  measurement  error  as 

a  reference. 


The  transition  between  SLM  and  a  constant  speed  turn  sees  a  reduction  in  error 
when  using  an  IMM  filter,  but  the  magnitude  is  still  significant.  The  combined  estimates 
distance  error  lies  directly  on  top  of  the  distance  error  for  model  two  in  Figure  9.  The 
estimates  take  approximately  thirty- five  samples  (0.35  seconds  with  a  0.01  sampling  rate) 
to  return  below  measurement  error  and  another  40  samples  (0.4  seconds)  before  returning 
to  steady-state.  The  transition  back  to  SLM  duplicates  the  results  seen  in  Figure  9.  A  later 
comparison  with  the  Kalman  filter  proves  these  transition  times  in  both  directions  to  be 
very  reasonable. 
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Figure  9.  Mean  distance  error  of  IMM  estimates  and  mean  measurement  error  during  a 
constant  speed  maneuver  initiated  at  sample  time  k  =  2000 


The  distance  errors  associated  with  a  constant  speed  turn  replicate  those  seen 
during  straight-line  motion  and  exhibit  only  half  of  a  meter  up-shift  towards 
measurement  error.  A  combination  of  transition  probabilities  (  pn  =  0.2  and 
p2X  =  0.4)  favoring  model  one,  IMM  based  prediction  methods,  and  high  process 
noise  (q2  =  4000)  during  the  turn  compared  to  SLM  (q  =  10)  take  away  from  model 
two  estimates  replicating  a  true  turn.  These  effects  also  drive  the  model  probability 
for  model  two  down.  The  result  is  seen  in  Figure  10  as  model  one  retains  the  lowest 
distance  errors. 
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Figure  10.  Constant  acceleration  mean  distance  error  of  IMM  estimates  with  mean 

measurement  error  as  a  reference. 


Another  useful  tool  in  assessing  IMM  performance  is  the  model  probability 
values  and  variations.  These  parameters  give  a  clear  indication  of  how  well  the  IMM 
reflects  reality.  The  probabilities  shown  in  Figure  11  are  averaged  over  100  samples  due 
to  the  vast  amount  of  model  probability  value  switching  that  occurred,  particularly  during 
the  turn.  Prior  to  the  turn,  model  one  probability  was  greater  but  not  drastically.  At  the 
onset  of  the  turn  (at  sample  2,000,  displayed  in  Figure  1 1  as  sample  20),  model  two  is 
assured.  However,  as  the  transition  to  steady-state  occurs,  the  IMM  becomes  less  certain 
that  a  turn  is  occurring.  The  transition  back  to  SLM  is  falsely  identified  as  an  assured 
target  turn,  when,  in  fact,  completion  of  the  turn  is  being  identified. 
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Figure  1 1 .  Model  probability  averaged  over  10  samples  during  IMM  mixing  with  the 
following  transitional  probabilities:  pu  =  0.8,/?12  =  0.2 ,p2l  =  ()AS,pu  =  0.55  . 

The  transitional  probabilities  used  in  the  simulation  described  above  were 
modified  based  on  the  inaccuracy  of  those  original  estimated  in  Equation  (40).  A  second 
simulation  is  presented  to  show  the  significance  of  the  adjustment.  The  only  parameters 
adjusted  from  the  above  simulation  were  pn  and  pn.  We  assumed  that  maneuvering 
would  be  minimal  and  estimated  the  values  accordingly  ( p2}  =  0.6  and  p22  =  0.4  ).  The 
simulation  resulted  in  higher  distance  errors  in  each  of  the  estimates  during  the  turn,  as 
seen  in  Figure  12.  The  combined  estimates  became  dominated  by  model  one.  The  SLM 
and  peak  model  transition  errors  were  unaffected,  but  an  analysis  of  the  model 
probabilities  indicated  a  major  problem  that  can  be  seen  in  Figure  13.  The  IMM  assumed 
the  target  remained  in  model  one  during  the  entire  simulation  with  the  exception  of  the 
transitions. 
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Figure  12.  Constant  speed  turn  mean  distance  error  of  IMM  estimates  with  mean 

measurement  error  as  a  reference. 
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Figure  13.  Model  probability  averaged  over  10  samples  during  IMM  mixing  with  the 
following  transitional  probabilities:  pn  =  0.8 ,pu  =  0.2, p21  =  0.6, p22  =  0.4. 
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Further  adjustments  to  the  transition  probabilities  proved  only  marginally 
effective.  The  values  used  in  the  first  simulation  were  adequate  to  track  the  missile 
through  a  turn.  We  assessed  the  IMM  performance  to  be  sufficient  and,  when  paired  with 
the  JPDAF,  capable  of  tracking  the  missile  during  split  target  tracking.  All  future 
simulations  will  ignore  the  model  one  and  model  two  estimates  and  be  evaluated  based 
on  combined  estimates  only. 

2.  Split  Target  Detection  and  Tracking 

True  split  target  tracks  and  model  parameters  were  adjusted  in  the  following 
simulation  to  display  split  track  generation  ability.  The  conditional  probabilities,  process 
and  sensor  noise,  sampling  frequency,  sample  length,  and  number  of  simulation  runs 
remain  unchanged.  The  generation  of  a  true  split  track  drives  model  three  to  generate 
measurements  for  the  split  and  attempt  track  formation.  The  transition  between  models 
three  and  four  is  the  focus,  and  a  parameter  t  represents  the  number  of  samples  model 
three  processes,  initially  set  to  five. 


Figure  14.  True  target  motion  focus  area  for  follow-on  figures  in  analysis  of  target  split 

simulations. 
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The  second  simulation  begins  in  SLM  with  target  splits  occurring  at  sample 
k  =100  and  at  sample  k  =800.  The  entire  simulation  of  1,000  samples  (sampling  rate 
remains  at  100  samples  per  second)  is  contained  within  Figure  14.  The  locations  in  which 
the  follow-on  figures  are  taken  are  also  displayed  in  Figure  14. 

The  convergence  of  the  filter  estimates  on  the  split  target  is  shown  in  Figure  15. 
Although  an  unlikely  decoy  deployment  scenario,  the  simulated  split  shows  the  ability  of 
the  filter  to  correct  to  a  drastic  change  of  original  estimates  (nearly  90  degrees  in  this 
case).  The  transition  to  the  red  line  corresponds  to  the  transition  to  model  four.  The  lower 
tracks  represent  the  missile  and  combined  estimates  and  are  where  model  three  estimates 
begin.  The  results  contained  in  Figure  16  show  an  enlargement  of  the  conditions  in 
Figure  15  at  the  split  junction.  Clear  divergence  of  the  model  three  estimates  occurs  after 
three  samples.  Model  four  takes  over  track  maintenance  and  is  indicated  by  the  solid  blue 
line  following  the  green  line  in  Figure  16. 


Figure  15.  Split  target  estimates  converging  on  split  target  from  SLM. 
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Figure  16.  Model  three  and  model  four  transitions  during  split  target  detection  from  SLM. 


The  simulation  continued  by  implementing  a  turn  from  sample  k  =  200  to 
k  =  900  and  generating  a  second  split  at  k  —  800.  The  same  color  scheme  as  above  is 
applied  to  the  true  and  estimated  tracks  for  the  split  a  from  constant  speed  maneuver.  The 
small  angle  of  separation  shown  in  Figure  17  allows  the  split  estimates  to  converge  more 
quickly  on  the  true  target  track.  On  the  other  hand,  the  small  divergence  of  the  split  track 
causes  the  model  three  estimates  to  not  noticeably  diverge  in  Figure  18  until  about  the 
fifth  sample. 
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Figure  17.  Split  target  estimates  converging  on  split  target  from  constant  acceleration. 


Figure  18.  Model  three  and  model  four  transitions  during  split  target  detection  from 

constant  acceleration. 
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The  stacked  vector  method  used  in  model  three  causes  a  linear  increase  in  mean 
distance  error  over  a  short  period  of  samples,  as  seen  in  Figure  19.  We  experienced  three 
major  factors  based  on  simulations  that  affect  the  magnitude  of  model  three  error:  the 
angle  of  separation  between  the  true  target  track  and  the  split  target  track,  the  missile 
state,  and  the  relative  speed  of  the  split  to  the  missile.  The  model  three  estimates  are 
passed  to  model  four,  which  receives  an  additional  measurement  and  performs  a  Kalman 
filter  update.  Hence,  the  maximum  mean  distance  error  seen  from  model  three  is  reduced 
in  one  estimate  to  the  initial  mean  distance  error  seen  in  Figure  19.  Both  split  tracks 
require  about  25  samples  to  match  measurement  error  and  another  25  samples  to  reach  a 
common  steady-state  error  value,  totaling  a  half  of  a  second  from  receipt  of  model  three 
estimates  to  nearly  one-quarter  of  measurement  error. 


Figure  19.  SLM  and  constant  acceleration  mean  distance  errors  during  target  splitting 

track  formation. 
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Figure  20.  SLM  and  constant  acceleration  mean  distance  errors  observed  during  split 

target  track  maintenance. 


The  combined  estimate  performance  is  reduced  as  a  result  of  the  model  three 
stacked  vector  method.  Although  difficult  to  ascertain  from  the  track  layout  during  the 
simulation,  the  mean  distance  errors  shown  in  Figure  20  show  a  20-meter  increase  for 
straight-line  motion  and  thirty-five  meter  increase  for  constant  acceleration.  The 
combined  estimates  require  approximately  60  samples  to  return  to  steady-state  following 
each  split  evolution.  The  steady-state  errors  are  consistent  on  either  side  of  the  respective 
splits;  however,  the  mean  distance  error  during  the  turn  is  raised  to  only  half  of  a  meter 
below  measurement  error  on  average.  The  model  one  and  model  two  transition  errors 
seen  in  Figure  9  are  duplicated  in  Figure  21. 
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Figure  2 1 .  Combined  estimate  mean  distance  error  versus  measurement  error  seen  during 

target  splitting. 


The  next  simulation  analyzes  the  effects  of  raising  the  number  of  measurements 
that  model  three  filters  before  transitioning  to  model  four.  For  a  split  target  track  that 
remains  significantly  close  to  the  true  target  track  (i.e.,  within  measurement  error),  further 
model  three  analyses  may  be  required  to  ensure  sufficient  track  separation  before  split 
target  confirmation.  The  recorded  common  past  history  maintained  by  the  stacked 
covariance  matrix  would  likely  be  more  accurate  than  relying  on  the  JPDA  with  no  cross¬ 
correlation  between  the  measurements. 

The  number  of  measurements  analyzed  by  model  three  t  was  raised  to  fifteen  with  all 
other  parameters  being  the  same  as  the  t  =  5  simulation.  Track  formation  in  Figure  22 
(analogous  in  position  to  Figure  16)  and  Figure  23  (analogous  in  position  to  Figure  18) 
are  improved  but  increasingly  lags  the  true  split  vector.  The  mean  distance  error 
continues  to  increase  linearly  as  seen  in  Figure  24.  The  time  required  for  the  split 
estimates  to  converge  on  the  split  target  were  minimally  effected  and  resulted  in  plots 
identical  to  Figure  15  and  Figure  27.  Therefore,  to  minimize  the  linearly  expanding  error 
and  still  be  able  to  generate  a  track,  a  value  of  t  =  2  will  be  used. 
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Figure  22.  Target  split  during  SLM  with  a  raised  number  of  model  three  samples  ( t  =  15 ). 


x  1 0 4 


Figure  23.  Target  split  during  constant  acceleration  with  a  raised  number  of  model  three 

samples  (  /  =  15). 
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Figure  24.  Model  three  mean  distance  error  versus  a  greater  number  of  samples  for  each 

target  split. 


C.  COMPARATIVE  ANALYSIS  WITH  A  STANDARD  KALMAN  FILTER 

1.  Kalman  Filter  Setup 

A  standard  Kalman  filter  was  created  and  run  in  parallel  with  the  IMMJPDA  filter 
for  split  target  comparison  purposes.  The  comparison  is  intended  to  evaluate  if  the 
additional  computing  of  the  IMMJPDA  results  in  increased  split  detection  performance 
over  a  basic  algorithm.  All  parameters  of  the  Kalman  filter  are  completely  independent 
from  the  IMM  and  JPDA  operations.  The  Kalman  filter  discriminates  between  the 
measurements  by  choosing  the  one  with  the  smallest  distance  from  the  predicted  track 
position.  Mean  distance  error  is  calculated  following  state  and  covariance  update  using 
the  norm  between  the  estimated  value  and  the  measurement.  The  plant  noise  was  adjusted 
to  allow  suboptimal  performance  during  SLM  in  exchange  for  reasonable  tracking 
capability  during  target  maneuver.  The  Kalman  can  be  adjusted  to  perform  better  in  either 
of  the  target  states  but  not  both.  The  steady-state  mean  distance  errors  can  be  seen  for 
SLM  in  Figure  25  and  for  a  constant  velocity  maneuver  in  Figure  27.  The  peaks  of  the 
distance  errors  in  Figure  26  were  matched  as  closely  as  possible  for  comparable 
maneuvering  tracking. 


45 


B 


8 

i 


(0 

o> 


14 

12 

10 


2 


50 


100 


Combined  Estimates 
-Kalman  Estimates 
"Measurement  Error 


150  200 

Sample  (k) 


250 


Figure  25.  Comparison  of  Kalman  and  IMMJPDA  filter  mean  distance  errors  during 

steady  state  SLM. 


Figure  26.  Comparison  of  Kalman  and  IMMJPDA  filter  mean  distance  error  peaks 
during  transition  from  SLM  to  constant  acceleration. 
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Figure  27.  Comparison  of  Kalman  and  IMMJPDA  filter  mean  distance  errors  during  a 

constant  acceleration  turn. 


2.  Scenario  Setup 

The  scenarios  designed  for  comparison  are  intended  to  replicate  decoy 
deployment  and  follow-on  ballistic  missile  evasion  tactics.  The  exo-atmospheric 
environment  allows  decoys  to  replicate  true  target  motion.  The  true  target  and  split  target 
vectors  are  presented  and  the  results  are  summarized  in  Table  2.  A  single  simulation 
consists  of  1,000  samples.  Each  scenario  is  used  to  generate  100  Monte  Carlo  simulation 
runs.  The  filters  are  evaluated  based  on  their  ability  to  maintain  true  target  track. 

Two  performance  criteria  are  used.  First,  if  the  final  mean  distance  error  value  of 
the  simulation  exceeds  a  threshold  value  (approximately  five  times  the  peak  values  seen 
in  Figure  26),  we  assume  the  decoy  has  tricked  the  filter.  The  percentage  of  success  out 
of  100  runs  is  annotated  in  Table  2.  Secondly,  a  counter  is  established  to  determine  when 
each  filter  chooses  split  target  measurements  over  true  target  measurements,  representing 
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the  filter  being  tricked  by  the  decoy  at  that  time  step.  The  measurements  associated  with 
the  filter  tracking  the  decoy  are  removed.  The  remaining  number  of  faulty  measurements 
associations  represents  isolated  instances  during  successful  target  tracking  where  decoy 
measurement  association  was  chosen.  Two  splits  occur  in  each  scenario.  Successful 
tracking  percentage  and  number  of  faulty  measurement  associations  are  recorded  per 
simulation  (1,000  samples)  for  each  split  separately.  In  general,  divergence  from  the  true 
target  is  about  20  degrees  for  the  first  split  and  by  less  than  10  degrees  for  the  second 
split.  The  distinguishing  features  between  the  scenarios  are  the  split  locations  and  the 
varying  target  state  before  and  after  the  split. 

Scenario  one  is  SLM  with  target  splits  at  k  =  200  and  k  =  600.  The  second  split 
is  intended  to  parallel  the  true  target  track.  The  layout  of  true  target  tracks  is  illustrated  in 
Figure  28. 


Figure  28.  Scenario  one  for  comparison  between  the  Kalman  and  IMMJPDA  filters. 
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The  second  scenario  contains  the  most  likely  method  of  decoy  deployment  and 
missile  maneuver.  The  split  occurs  at  k  =200,  followed  by  a  turn  from  k  =  205  to 
k—  500.  The  second  decoy  is  deployed  at  k  =  600  and  followed  by  a  turn  from  k  —  605 
to  k  =  1000 .  The  true  target  and  split  vectors  are  shown  in  Figure  29. 


Figure  29.  Scenario  two  for  comparison  between  the  Kalman  and  IMMJPDA  filters. 


Although  a  constant  missile  maneuver  is  detrimental  to  fuel  consumption  and 
increases  vulnerability  to  tracking  systems,  the  third  scenario  is  included  for  a  complete 
analysis  of  each  target  state.  The  missile  deploys  decoys  and  maintains  constant  velocity 
throughout  the  simulation  while  maneuvering,  as  shown  in  Figure  30. 
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Figure  30.  Scenario  three  for  comparison  between  the  Kalman  and  IMMJPDA  filters. 


The  fourth  scenario  tests  the  feasibility  of  tracking  a  decoy  during  the  transition 
from  SLM  to  constant  speed  acceleration.  The  first  split  occurs  10  samples  following  the 
commencement  of  a  turn.  The  ten-sample  delay  (instead  of  a  25 -sample  delay)  was 
chosen  to  prevent  an  advantage  to  either  filter  and  was  chosen  based  on  the  results  of 
Figure  26.  The  turn  is  maintained  for  200  samples  before  returning  to  SLM.  The  second 
split  is  deployed  under  the  same  conditions  but  with  significantly  less  divergence  from 
the  missile.  The  true  target  vectors  can  be  seen  in  Figure  31. 
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Figure  3 1 .  Scenario  four  for  comparison  between  the  Kalman  and  IMMJPDA  filters. 

Similar  to  the  previous  conditions,  the  fifth  scenario  analyzes  the  ability  to  track  a 
decoy  through  a  transition  from  constant  velocity  acceleration  to  SLM.  Distance  error 
spikes  occur  during  model  transitions  despite  the  direction  as  seen  in  Figure  21.  The 
missile  begins  in  a  turn  and  transitions  to  straight-line  motion  followed  by  a  split  10 
samples  later.  The  second  split  is  under  similar  conditions,  and  both  split  locations  can  be 
seen  in  Figure  32. 
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Figure  32.  Scenario  five  for  comparison  between  the  Kalman  and  IMMJPDA  filters. 


The  final  scenario  involves  the  missile  maintaining  SLM  throughout  and 
deploying  both  splits  within  five  samples  of  each  other.  The  first  split  is  deployed  from 
the  left  side  of  the  missile  at  £  =  500,  and  the  second  is  ejected  at  k  =  505  from  the  right 
side  of  the  missile.  The  vicinity  of  both  splits  to  the  true  target  is  identical  to  the  second 
split  in  scenario  one.  No  visual  plot  is  included  for  this  scenario,  but  the  results  are 
included  in  Table  2.  The  results  for  this  scenario  are  expressed  as  a  single  event  and  are 
not  differentiated  between  the  two  splits. 

3.  Results 

The  IMMJPDA  filter  displayed  improved  missile  tracking  performance  and  decoy 
avoidance  over  the  tuned  standard  Kalman  filter  during  the  specified  scenarios.  The 
results  between  the  two  splits  of  the  first  scenario  show  the  benefits  of  the  JPDA  filter 
over  a  simple  nearest  neighbor  algorithm. 

The  second  and  fourth  scenarios  validate  the  additional  computations  required  of 

the  IMMJPDAF  over  the  standard  Kalman  filter.  We  expected  a  significant  drop  in  the 
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ability  to  maintain  true  target  track  given  the  large  distance  error  during  model  switching 
in  comparison  to  sensor  error.  The  successful  tracking  of  the  true  target  by  the 
IMMJPDAF  in  scenario  two  and  scenario  four  discredit  the  concerns  about  IMMJPDAF 
degradation  between  model  transitions.  The  number  of  faulty  measurement  associations 
in  both  scenarios  highlights  a  need  for  model  improvements.  However,  these  results 
prove  that  the  combination  of  the  IMM  and  the  JPDA  can  compensate  for  simple  target 
models. 

The  ability  to  track  a  decoy  deployment  through  a  turn  is  proven  in  scenario  three. 
A  repeating  trend  throughout  each  simulation  is  the  effect  of  the  non-acting  model 
estimate  upon  the  existing  model.  The  increased  number  of  faulty  samples  over  the 
standard  Kalman  filter  is  explained  by  the  model  one  estimates  that  tend  to  associate 
more  with  the  split  target  in  this  scenario.  The  flawless  success  percentage  proves 
adequate  model  transition  probabilities  to  allow  the  filter  to  self-correct  with  the  model 
two  estimates. 

Despite  scenario  five  being  a  non-traditional  decoy  deployment  strategy,  the 
IMMJPDAF  performed  as  expected  and  similarly  to  scenarios  two  and  four.  The  results 
of  these  conditions  show  consistent  performance  through  model  transitions. 

The  results  of  the  final  scenario  suggest  the  IMMJPDAF  may  be  suitable  for 
ICBM  tracking  in  a  high  decoy  density  scenario.  The  decoys  maintained  essentially  a 
constant  track  (each  velocity  component  was  offset  by  merely  0.01  meters  per  second  in 
each  direction  for  both  decoys)  to  the  missile,  yet  tracking  accuracy  was  maintained.  The 
combination  of  target  track  gates  and  hypothesis  scoring  yielded  significant 
improvements,  especially  when  compared  to  the  standard  Kalman  filter. 
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Table  2.  Performance  summary  of  the  Kalman  filter  and  the  IMMJPDAF  under 
simulated  ICBM  decoy  deployment  and  maneuver  scenarios.  Refer  to 
section  two  (Scenario  Setup)  for  an  explanation  of  table  organization. 


Scenario 

Split  Number 

Kalman  Filter 

IMMJ] 

PDA  Filter 

Success 

(%) 

Faulty 

Measurement 

Associations 

Per 

Successful 

Simulation 

Success 

(%) 

Faulty 

Measurement 

Associations 

Per 

Successful 

Simulation 

1 

1 

100 

0.25 

100 

0.14 

2 

71 

1.3 

100 

0.35 

2 

1 

0 

100 

254.6 

2 

0 

100 

22.1 

3 

1 

100 

0.04 

100 

0.21 

2 

100 

0 

100 

0.16 

4 

1 

0 

100 

23.1 

2 

0 

100 

147.8 

5 

1 

0 

100 

122.1 

2 

0 

100 

176.4 

6 

1 

59 

0.1 

100 

0.71 

2 

54 


V.  CONCLUSIONS 


A.  SUMMARY  OF  RESULTS 
1.  IMM 

With  regard  to  ICBM  tracking,  the  IMM  proved  to  be  capable  during  both  steady- 
state  and  maneuvering  conditions.  For  the  given  conditions,  the  estimates  resulting  from 
the  IMM  mixing  process  received  a  reduction  from  measurement  error  of  50  percent  on 
average.  The  degradation  of  the  individual  model  estimates  was  minimal.  A  significant 
amount  of  iterations  was  required  to  find  optimal  values  for  transition  probabilities  ptj. 

Even  with  the  determined  values  for  ptj,  the  model  probability  «. transitions,  especially 

during  a  constant  velocity  maneuver,  were  frequent  and  often  occurred  between  each 
sample.  For  targets  with  the  likelihood  of  transitioning  between  established  models  often, 
this  method  is  appropriate.  However,  for  the  purposes  of  an  ICBM  that  prefers  deploy 
decoys  to  conducting  maneuvers,  the  necessity  to  account  for  continuous  contact 
maneuver  is  reduced.  In  addition,  the  distance  errors  seen  during  model  transitions  are 
concerning  due  to  the  likelihood  of  decoy  deployment  during  those  periods.  The 
IMMJPDA  filter  performed  well  during  these  transitions  for  the  scenarios  tested  but 
could  become  a  concern  under  higher  split  target  densities. 

With  regard  to  split  target  tracking,  the  IMM  provided  a  convenient  framework  in 
which  to  generate  a  split  target  and  allowed  for  a  seamless  transition  for  follow-on  track 
maintenance.  We  introduced  a  non- traditional  IMM  approach  where  the  two  models 
regarding  the  split  target  are  not  subject  the  standard  IMM  processes.  Initially  stacking 
the  combined  estimates  to  generate  the  split  track  allowed  the  common  history  of  the 
target  and  decoy  to  be  exploited.  The  decoy  track  generation  caused  an  increase  in 
distance  error  of  the  combined  estimates  of  over  two  and  a  half  times  measurement  error. 
The  decoupling  of  model  three  and  model  four  with  the  target  models  prevented  further 
degradation  of  the  combined  estimates.  Additionally,  the  linearly  increasing  model  three 
errors  showed  that  decoupling  of  the  split  target  and  true  target  sooner  was  best. 
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Prolonged  tracking  of  the  split  target  was  detrimental  to  both  the  combined  estimates  and 
the  split  track  maintenance. 

2.  JPDA 

The  JPDA  was  able  to  discriminate  between  true  target  and  decoy  measurements 
that  fell  within  established  track  gates  with  minimal  faulty  associations.  Although 
dependent  on  the  quality  of  estimates  provided  and  difficult  to  evaluate  on  its  own,  the 
JPDA’s  comparison  to  the  nearest  neighbor  algorithm  implemented  by  the  Kalman  filter 
proved  very  effective.  The  generation  of  hypothesis  is  simplistic  for  the  scenarios 
presented,  but  computational  requirements  drastically  increase  with  additional 
measurements. 

B.  RECOMMENDATIONS  FOR  FUTURE  WORK 

1.  Probability  of  Detection  and  False  Returns 

The  assumption  that  the  probability  of  detection  PD  is  one  requires  a  sensor  that 
produces  a  valid  measurement  for  each  target,  including  all  decoys,  at  every  measurement 
time.  The  assumption  on  no  false  alarms  (e.g.,  clutter  or  spurious  measurements)  reduces 
the  number  of  hypotheses  that  the  JPDA  filter  must  consider.  Measurements  may  be 
included  that  do  not  belong  to  any  track,  or  an  established  track  may  not  receive  an 
associated  measurement  at  any  given  time  interval.  The  expansion  of  the  JPDA  filter 
equations  to  include  these  parameters  is  discussed  in  [3],  This  would  result  in  a  more 
complicated  but  potentially  more  flexible  JPDA  filter  and  is  worth  future  exploration. 

2.  Model  Modifications 

The  established  models  to  describe  the  ICBM  dynamics  do  not  fully  describe  a 
body  traveling  in  exo-atmospheric  conditions.  An  inclusion  of  the  gravitational  forces  on 
the  missile  on  the  z  -dimension,  rotating  Earth  effects  [2],  and  the  associated  filter 
performance  should  be  explored. 

As  decoys  become  more  sophisticated,  the  expansion  of  a  decoy  model  beyond 

mimicking  the  trajectory  of  the  missile  may  be  necessary.  Application  of  a  separate  IMM 
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process,  complete  with  its  own  set  of  models,  for  a  split  target  may  be  appropriate  for 
splits  determined  to  be  a  reentry  vehicles. 

3.  Algorithm  Flexibility 

Given  the  existing  MATLAB  code  as  a  guide,  expansion  to  include  any  number 
of  splits  would  be  essential  for  further  simulation  testing.  As  is,  the  algorithms  can 
support  zero,  one,  or  two  splits.  Any  additional  splits  written  in  the  same  manner  will 
require  an  increasing  amount  of  code  for  each  one.  The  additional  hypothesis  required  for 
the  JPDA  and  the  gate-related  calculations  contribute  to  the  majority  of  the  required  code 
expansion. 
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APPENDIX, 


MATLAB  CODE 


A.  MASTER  SCRIPT  FILE  INCLUDING  ALL  INITIAZATIONS,  ORDER  OF 
FUNCTION  FILE  CALLS,  AND  GENERATION  OF  ALL  PLOTS 

close  all; 
clear  all; 
clc; 

delta=0 . 01 ; 
nsamples=1000  ; 
nloops=100 ; 

%Turn  start  and  stop  times 
t_startl=300  ; 
t_stopl=  700; 
t_start2=1400 ; 
t_stop2=l 600 ; 

%Split  locations 
n= [1200; 1503] ; 

mdl= [ ] ; 
md2=[]  ; 
mdo=  [  ] ; 
md3=  [  ] ; 
md4a=  [  ] ; 
md4b=  [  ] ; 
mdk=  [  ]  ; 

for  ii=l:nloops 

%Measurement  matrix 
H=  [  1  0  0  0  0  0 ; 

0  0  1  0  0  0; 

0  0  0  0  1  0]; 

%Sensor  Error 
se=5  ; 

sigmax  =  se; 
sigmay  =  se; 
sigmaz  =  se; 

R=diag ( [ sigmax A 2 ; sigmayA2 ; sigmaz A 2 ] ) ; 

%Target  state  transition  matrix 
FI  =  [1  delta  0  0  00; 

0  1  0  0  0  0; 
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0  0  1  delta  00; 

0  0  0  1  0  0; 

0  0  0  0  1  0; 

0  0  0  0  0  1]; 


%Initialize  the  target  vectors 


[ 


0; 

2500; 

0; 

2300; 
10000; 
0  ]  ; 


PI  =  10 A4 . *  [100000; 

0  1  0  0  0  0; 

0  0  1  0  0  0; 

0  0  0  1  0  0; 

0  0  0  0  1  0; 
000001]; 


%x  position 
%x  velocity  in  m/ s 
%y  position 
%y  velocity  in  m/ s 
%z  position 

%z  velocity  in  m/s 

%Covariance  matrix 


ql=l 0 ; 

Ql  =  ql 

0; 

0; 

0; 

0; 


[ (deltaA3 ) /3  (deltaA2) /2 


(deltaA2 ) /2 
0 
0 
0 

(deltaA3)/3  (deltaA2)/2; 
0 


(deltaA2 ) / 2 


Ql- 


0 

ql2 

0 

0 

0 

0 


delta] 

0 

0 

0 

0 

0 

0 


ql2  0 
0  0 
0  0 


delta 

0 

0 

0 

0 


0; 

0; 

0; 

0; 

0; 

ql2; 


0  0  0 

0  0  0 

(deltaA3 ) / 3  (deltaA2)/2  0 

(deltaA2 ) / 2  delta  0 

0  0 

0  0 


Qk=200*Ql; 


spd=  3397.06; 
g=-20 ; 


%in  m/ s 
%  t  of  g's 


omega=g*  9 . 80  665/spd; 

F2=  [1  sin (omega*delta) /omega 


%turn  rate  in  radians/sec 
0  ( 1-cos (omega*delta) ) /omega  0 
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0; 


0 


0; 

0 

cos (omega*delta) 

0; 

0 

(1-cos (omega*delta) ) /omega 

0; 

0 

sin (omega*delta) 

0; 

0 

0 

1]  ; 

0 

0 

P2 

=  10A 

4  .  * 

[1  o 

0 

0 

0 

0; 

0  1 

0 

0 

0 

0; 

0  0 

1 

0 

0 

0; 

0  0 

0 

1 

0 

0; 

0  0 

0 

0 

1 

0; 

0  0 

0 

0 

0 

1]  ; 

q2- 

o 

o 

o 

■v 

Q2 

=  q2 

.*  [ (deltaA 

3)  / 3 

(deltaA2 ) /2 

0; 

(delta 

A2 )  / 2 

delta 

0; 

0 

0 

0; 

0 

0 

0; 

0 

0 

0; 

0 

0 

1]  ; 

F3 

=  [1 

delta 

0 

0 

0 

1 

0 

0 

0 

0 

1 

delta 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

Q3 

=  ql 

.*  [ (deltaA 

3)  / 3 

(deltaA2 ) /2 

0; 

(delta 

A2 )  / 2 

delta 

0; 

0 

0 

0; 

0 

0 

0/ 

0 

0 

0/ 

0 

0 

0  -sin (omega*delta) 

1  sin (omega*delta) /omega  0 
0  cos (omega*delta)  0 
0  0  1 
0  0  0 


%Covariance  vector 


0  0  0 
0  0  0 

(deltaA3 ) / 3  (deltaA2)/2  0 

(deltaA2)/2  delta  0 

0  0  1 

0  0  0 


0  0; 

0  0; 

0  0; 

0  0; 

1  0; 

0  1]  / 


0  0  0 
0  0  0 

(deltaA3 ) / 3  (deltaA2)/2  0 

(deltaA2)/2  delta  0 

0  0  1 

0  0  0 
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P3  =  10A4.*  [100000;  %Covariance  vector 

0  1  0  0  0  0; 

0  0  1  0  0  0; 

0  0  0  1  0  0; 

0  0  0  0  1  0; 

000001]; 


K1  =  P1*H* *inv (H*P1*H'  +  R) ; 
K3  =  P3*H' *inv (H*P3*H'  +  R) ; 
Ks  =  [K1  zeros (6,3) ; 

zeros (6,3)  K3 ] ; 


0  0; 

0  0; 

0  0; 

0  0; 

1  0; 

0  1]; 


F4  = 


[1 

0 

0 

0 

0 

0 


delta 

1 

0 

0 

0 

0 


0 

0 

delta 

1 

0 

0 


Q4 

0; 

0; 


ql  .*  [  (deltaA3 ) / 3 

(deltaA2 ) /2 
0 


0; 


0 


0; 


0 


0; 


0 


1]  ; 


(deltaA2 ) /2 


delta 

0 

0 

0 

0 


0  0 
0  0 


(deltaA3 ) /3  (deltaA2) /2 
(deltaA2)/2  delta 
0  0 

0  0 


[1 

0 

0 

0 

0 

0; 

0 

1 

0 

0 

0 

0; 

0 

0 

1 

0 

0 

0; 

0 

0 

0 

1 

0 

0; 

0 

0 

0 

0 

1 

0; 

0 

0 

0 

0 

0 

1]  ; 

%Covariance  vector 


xlh=x; 

x2h=x; 

xk=x; 

Pk=Pl; 

x4 a= [ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ]  ; 
x4b= [ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ] ; 
tcorrect=0 ; 


xo=x; 


0 


0 

0 

0 

1 

0 
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Po=Pl; 


%Intial  Likelihood  for  each  state 

ul  =  0.7;  %More  likely  to  remain  in  straight  line  motion 
u2  =  0.3;  %More  likely  to  transition  out  of  turn  and  int 
motion 

%Probability  of 
pll  =  0.75; 
pl2  =  0.25; 
p2 1  =  0.45; 
p22  =  0.55; 

%Model  1  and  Model  2  mixing 
clb  =  pll*ul  +  p21*u2; 
c2b  =  pl2*ul  +  p22*u2; 

%Matrix  Initialization 
xlhm= [ ] ; 
znl= [  j ; 
eml= [ ] ; 
x2hm= [ ] ; 
zn2  = [ ] ; 
em2= [ ] ; 
x3hm= [ ] ; 
zn3= [ ] ; 
em3= [ ] ; 
xsh=  [  ]  ; 
zn4- [ j ; 
x4hm= [ ] ; 
em4a= [ ] ; 
em4b= [ ] ; 
xm= [ ] ; 
x3m=  [  ]  ; 
x4am= [ ] ; 
x4bm= [ ] ; 
xom=  [  ]  ; 
znk= [ ] ; 
xkm= [ ] ; 
dlm= [ ] ; 
d2m= [ ] ; 
dkm= [ ] ; 
dom= [ ] ; 
d3m= [ ] ; 
d4am= [ ] ; 
d4bm= [ ] ; 
znkm= [ ] ; 
znc=zeros  (6,1)  ; 
mitrick= [ ] ; 
mktrick= [ ] ; 


changing  state 

%Probability  of  changing  state 


than  to  turn 
stright  line 
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%samples  in  model  3  before  transition  to  4 
t=l ; 


for  k= 1 : n  s amp 1 e  s 


%Measurement  matrix 

H= [1  0  0  0  0  0 ; 

0  0  1  0  0  0; 

000010]; 

%Sensor  Error 
sigmax  =  se; 
sigmay  =  se; 
sigmaz  =  se; 

R=diag ( [ sigmax A 2 ; sigmayA2 ; sigmaz  A 2 ] )  ; 

%Measuring  Target  Position 

ztrue  =  H*x;  %True  position 

%Add  sensor  noise  to  the  measurements  and  calculate  measurement 

error 

w=randn ( 3 , 1 ). * [ sigmax;  sigmay;  sigmaz]; 

znu  =  ztrue  +  w; 

em  =  sqrt(w'*w); 

znl=[znl,  znu] ; 

eml= [eml , em] ; 


%Target  motion  step  for  split  target 


%Generate  the  first  split  vector 
if  k==n ( 1 ) 

x3  =  x  +  [0;  0 . 2  *  x  ( 2  )  ;  0;  +0.5*x(4);  0;  0]; 

x3h=xo ; 

xsh= [xo ; x3h] ; 

x3hm= [x3hmA  x3h] ; 

x3m= [x3m, x3 ] ; 

Ko  =  Po*H ' ^ inv (H^Po^H '  +  R) ; 

K3  =  P3*H ' *inv (H^P3^H '  +  R) ; 

Ks  =  [Ko  zeros(6,3); 

zeros (6, 3)  K3 ] ; 
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end 


Ps  =  [Po 
Po 


Po; 

Po]  / 


if  k==n (2 ) 

x3  =  x  +  [0;  0 . 1 * x  ( 2 )  ;  0;  -0.1*x(4); 

x3h=xo ; 

xsh= [xo ; x3h] ; 

x3hm= [x3hm,  x3h] ; 

x3m= [x3m, x3 ] ; 


o 

> 

* 

[1 

0 

0  0  0 

0; 

0 

1 

0  0  0 

0; 

0 

0 

10  0 

0; 

0 

0 

0  10 

0; 

0 

0 

0  0  1 

0; 

0 

0 

0  0  0 

1] 

r 

Ko  = 

Po*H  ' 

*inv (H*Po*H ' 

+ 

R)  ; 

K3  = 

P3*H ' 

*inv (H*P3*H ' 

+ 

R)  ; 

Ks  = 

[Ko 

zeros  ( 6 

,3) 

zeros (6,3) 

K3  ] 

Ps  = 

[Po 

Po; 

Po 

Po] 

r 

0/  0]/ 


%Covariance  vector 


end 


%Generate  true  split  motion  for  x3 
if  k  >  n(l)  &&  k<=n(l)+t 
x3=F3*x3 ; 
x3m= [x3m, x3 ] ; 

end 

if  k  >  n(2)  &&  k<=n(2)+t 

x3=F3*x3 ; 
x3m= [x3m, x3 ] ; 


end 


%Generate  true  target  position  for  x4,  two  splits 
if  k==n ( 1 ) +t 
x4a=x3 / 
x4= [x4a; x4b] ; 
x4am= [x4am,x4a] ; 


elseif  k>=n(l)+t 
x4a=F4,lrx4a; 
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x4= [x4a; x4b] ; 
x4am=[x4am,x4a] ; 

end 

if  k==n ( 2 ) +t 
x4b=x3 ; 
x4  = [x4a; x4b] ; 
x4bm= [x4bm, x4b] ; 

elseif  k>n ( 2 ) +t 
x4b=F4*x4b; 
x4  = [x4a; x4b] ; 
x4bm= [x4bm, x4b] ; 

end 


if  k>n ( 1 ) +t 

[  x4h, P4 ,  znu , znc, d4a , d4b , itrick  ]  =  modelfour( 
delta, x4h,x4,P4,Q4,xlh,H,R,se,  znc, znu, PI , P3  )  ; 
zn4  = [ zn4 ,  znc] ; 
x4hm= [x4hm,  x4h] ; 
d4am= [d4am,  d4a] ; 
d4bm= [d4bm,  d4b] ; 
mitrick=[ mi trick, itrick] ; 


else 


x4=zeros (12,1) ; 
x4h=zeros (12, 1) ; 

[  x4h, P4 , znu, znc, d4a, d4b, itrick  ]  =  modelfour( 
delta, x4h,x4,P4,Q4,xlh,H,R,se, znc,  znu,  PI ,  P3  )  / 
zn4  = [ zn4 , znc] ; 
x4hm= [x4hm, x4h] ; 
d4am= [d4am, d4a] ; 
d4bm= [d4bm, d4b] ; 
mitrick= [mi trick, itrick] ; 


end 

if  k  >n(l)  &&  k<=n(l)+t 

[PI, P3, Ps,Ks, zns,xlh,x3h,dl,d3, itrick]  = 
model three (delta, x,x3,xsh,Q3,Ps,H,R,F3,Ks, PI, se,  znu) / 
mitrick= [mi trick, itrick] ; 
zn3=[zn3,zns] ; 
x3hm= [x3hm, x3h] / 
dlm=  [dim,  dl  ]  ; 
d3m=  [d3m,  d3  ]  ; 
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[x2h,P2,d2]  =  modeltwo (x2h, x, P2 ,  znu, H, R) ; 
d2m= [d2m,  d2 ] ; 

znk= [ znu ; zns ; zeros (3,1) ] ; 
z  n  km= [znkm, znk] ; 

[xk, Pk, dk, ktrick]  =  kalman (  delta, xk, x, Pk, znk, Qk, se ) ; 
mktrick= [mktrick, ktrick] ; 
xkm= [xkm, xk] ; 
dkm= [dkm, dk] ; 


[xOlh, x02h, P01 , P02 , ul , u2 , xo, Po, do, clb, c2b] =modelprob (xlh,x2h,Pl,P2,ul,u 
2 ,  znu,  H, R, x, clb, c2b)  ; 

xom= [xom, xo ] ; 
dom= [dom, do] / 


[xlh, x2h, PI , P2 ] =prediction (xOlh, x02h, PO 1 , P02 , FI , F2 , Q1 , Q2 , ul , u2 ) ; 
xlhm= [xlhm, xlh] ; 
x2hm= [x2hm, x2h] ; 

elseif  k  >n(2)  &&  k<=n(2)+t 

[PI, P3, Ps,Ks, zns,xlh,x3h,dl,d3, itrick]  = 
model three (delta, x,x3,xsh,Q3,Ps,H,R,F3,Ks, PI, se, znu) / 
mitrick= [mi trick, itrick] ; 

zn3=[zn3, zns] / 
x3hm= [x3hm, x3h] ; 
dlm= [dim, dl ] ; 
d3m=  [d3m,  d3  ]  / 

[x2h,P2,d2]  =  modeltwo (x2h, x, P2 , znu, H, R) ; 
d2m= [d2m, d2 ] ; 

znk= [ znu; znc] ; 
z  n  km= [ z  n  km , znk] ; 

[xk, Pk, dk, ktrick]  =  kalman (  delta, xk, x, Pk, znk, Qk, se ) ; 
mktrick= [mktrick, ktrick] ; 
xkm= [xkm, xk] ; 
dkm= [dkm, dk] ; 


[xOlh, x02h, P01 , P02 , ul , u2 , xo, Po, do, clb, c2b] =modelprob (xlh,x2h,Pl,P2,ul,u 
2 ,  znu,  H, R, x, clb, c2b)  ; 
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xom= [xom, xo ] ; 
dom= [dom, do] ; 


[xlh,  x2h, PI , P2 ] =prediction (xOlh, x02h, PO 1 , P02 , FI , F2 , Q1 , Q2 , ul ,  u2  ) ; 
xlhm= [xlhm, xlh] ; 
x2hm= [x2hm, x2h] ; 

else 

[  xlh,Pl,dl,R]  =  modelone (  xlh, x, PI , H, R, znu) ; 
dlm= [dim, dl ] ; 

[x2h,P2,d2]  =  modeltwo (x2h, x, P2 , znu, H, R) / 
d2m= [d2m, d2 ] ; 

if  k<n ( 1 ) 

znk= [ znu; zeros (6,1) ] ; 

else 

znk= [ znu; znc] ; 

end 

z  n  km=  [  z  n  km ,  z  n  k  ]  ; 

[xk, Pk, dk, ktrick]  =  kalman (  delta, xk, x, Pk, znk, Qk, se ) ; 
mktrick= [mktrick, ktrick] ; 
xkm= [xkm, xk] ; 
dkm= [dkm, dk] ; 


[xOlh, x02h, P01 , P02 , ul , u2 , xo, Po, do, clb, c2b] =modelprob (xlh,x2h,Pl,P2,ul,u 
2 ,  znu,  H, R, x, clb, c2b)  ; 

xom= [xom, xo ] ; 
dom= [dom, do] ; 


[xlh, x2h, PI , P2 ] =prediction (xOlh, x02h, PO 1 , P02 , FI , F2 , Q1 , Q2 , ul , u2 ) ; 
xlhm= [xlhm, xlh] ; 
x2hm= [x2hm, x2h] ; 


end 


if  k==n ( 1 ) +t 


x4  ( 1 : 6, 1 ) =x3 ; 
x4  (7  : 12,l)=x4b; 

%x4h ( 1 : 6 , 1 ) =x3h; 
%x4h(7:12,l)=zeros(6,l) ; 
x4h= [x3h; zeros (6,1) ] ; 
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x4hm= [x4hm, x4h] ; 
znc=zeros (6,1) ; 

elseif  k==n(2)+t 

x4 (7 : 12, 1) =x3; 
x4h (7 : 12, 1) =x3h; 
x4hm= [x4hm, x4h] ; 


end 

%Target  Motion  Step  for  Straight  line  motion 
x=Fl*x; 

%Target  motion  step  for  constant  acceleration 
if  k>t_startl  &&  k<t_stopl 
x=F2  *x; 

elseif  k>t_start2  &&  k<t_stop2 
x=F2  *x; 

end 

%Collect  true  target  vectors 
xm= [ xm , x ] / 


end 


if  dkm ( size (dkm, 2 ) ) >60 0 
ttk=l ; 

else 

ttk=0 ; 

end 

if  dom (size (dom, 2 ) ) >600 
tto-1 ; 

else 

tto=0 ; 

end 


Keep  a  running  total  of  distance  errors 
f  ii==l 

mdl=dlm; 

md2=d2m; 

mdo=dom; 

md3=d3m; 

md4a=d4am; 

md4b=d4bm; 

mdk=dkm; 
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meml=eml ; 

txom=xom; 

tznl=znl; 

txkm=xkm; 

tzn3=zn3 ; 

tx3hm=x3hm; 

tzn4=zn4; 

tx4hm=x4hm; 

txkm=xkm; 

tottk=ttk; 

totto=tto ; 

sitrick=sum (mitrick) ; 
sktrick=sum (mktrick) ; 

else 

mdl=mdl+dlm; 

md2=md2+d2m; 

mdo=mdo+dom; 

md3=md3+d3m/ 

md4a=md4a+d4am; 

md4b=md4b+d4bm; 

mdk=mdk+dkm; 

meml=meml+eml ; 

txom=txom+xom; 

tznl— tznl+znl ; 

txkm=txkm+xkm; 

tzn3=tzn3+zn3; 

tx3hm=tx3hm+x3hm; 

tzn4=tzn4+zn4; 

tx4hm=tx4hm+x4hm; 

txkm=txkm+xkm; 

tottk=tottk+ttk 

sktrick=sktrick+sum (mktrick) 

totto=totto+tto 

sitrick=sitrick+sum (mi trick) 


end 


end 

%Calculate  average  distance  errors 

mdl=mdl/ nloops ; 

md2=md2 / nloops ; 

mdo=mdo/ nloops ; 

md3=md3/ nloops ; 

md4a=md4a/ nloops ; 

md4b=md4b/nloops / 

mdk=mdk/nloops ; 

meml=meml /nloops / 

mxom=txom/ nloops / 
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mznl=tznl/nloops ; 
mxkm=txkm/ nloops / 
mzn3=t zn3/nloops / 
mx3hm=tx3hm/ nloops ; 
mzn4=t zn4 /nloops / 
mx4hm=tx4hm/ nloops ; 
mxkm=txkm/ nloops; 
titrick=sitrick/nloops; 
tktrick=sktrick/nloops; 


figure ; 
hold  all; 
plot (mdo, ' g ' ) ; 
plot (mdk, '  r  '  )  ; 
plot (meml , '  b  '  )  ; 

title (' Combined (Model  1/2)  vs  Kalman  Estimate  Mean  Error'); 
set  (gcf , 'color' ,  '  w '  ) ; 

legend  ('Combined  Estimates ',' Kalman  Estimates ' ,  ' Measurement  Error'); 
grid  on; 

set (  f indob j (gca,  type ' ,  ' line ') ,  'LineWidth',  1.5); 
xlabel ( ' Sample  ( k)  ' ,  ' font size '  ,  18 ,  ' f ontweight ' ,  ' bold ' ) ; 
ylabel('Mean  Distance  Error 

(meters ) ' ,  ' font size ' , 18, ' f ontweight ' ,  ' bold ' ) ; 

set  (  findobj  (gca,  type ',' line ')  ,  'LineWidth',  1.5); 


figure ; 

title ( ' IMM  and  JDPA' ) 
hold  all; 

plot (xm ( 1 , : ) , xm (3 , : ) , ' b ' ) ; 
plot (xkm ( 1 , : ) , xkm ( 3 , : ) , ' r ' ) ; 
plot (xom ( 1 , : ) , xom ( 3 , : ) , ' m ' ) ; 

First  split 

plot(zn3(l,l:size(zn3,2)/2) , zn3(2,l*size(zn3,2)/2) , 'g') ; 
plot (x3m ( 1 , 1 : size (x3m, 2 ) / 2 ) , x3m ( 3 , 1 : size (x3m, 2 ) /2 ) , 'g') ; 
plot (x3hm ( 1 , 1 t size (x3hm, 2) / 2) ,  x3hm ( 3 , 1 : size (x3hm, 2) / 2)  ,  ' r ' ) ; 
plot(mzn4 ( 1 , 1 : size (mzn4 , 2 )  ) ,zn4  (2,l:size(zn4,2) ) ,  ' r '  )  ; 
plot (x4am ( 1 , 1 i size (x4am, 2 ) ) , x4am (3 , 1 : size (x4am, 2 ) ) , ' g ' ) ; 
plot (x4hm ( 1 , 1 t size (x4hm, 2 ) ) , x4hm (3 , 1 : size (x4hm, 2 ) ) , ' r ' ) ; 

Second  split 

plot (zn3 (1, size (zn3, 2) /2+l:size (zn3, 2) ) , zn3 (2, size (zn3, 2) /2+l:size (zn3, 
2) ) , 'g') ; 

plot (x3m ( 1 , size (x3m, 2)/2+l:size (x3m, 2 ) ) , x3m ( 3 , size (x3m, 2)/2+l:size (x3m, 
2) ) , 'g' ) 

plot (x3hm ( 1 , size (x3hm, 2)/2+l:size (x3hm, 2) ) , x3hm ( 3 , size (x3hm, 2)/2+l:size 
(x3hm, 2 ) ) , 'r') ; 

plot(mzn4 ( 4 , n (2 ) -n ( 1 ) +1 : size (mzn4 , 2 ) ) ,zn4 (5, n (2 ) - 
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n(l)+l:size(mzn4,2) ) , ’ r ’ ) ; 

plot (x4bm ( 1 , 1 : size (x4bm, 2 )  )  , x4bm (3 , 1 : size (x4bm, 2 ) ) ,  ' g ’ ) ; 
plot  (x4hm (7  ,  n (2) -n (1) +1 : size (x4hm, 2 ) ) , x4hm (  9,  n  (2  )  - 
n(l)+l:size (x4hm, 2 ) ) ,  ' r ' ) ; 

set  (gcf , 'color' , '  w '  ) ; 

legend  (’True  Missile  Track' , ' True  Split  Track ',' Split  Estimates'); 
grid  on; 

set (  f indob j (gca, ' type ',' line ') ,  'LineWidth',  1.5); 
xlabel ( ' X  (meters) ' , ' font size ' , 18, ' f ontweight ' , ' bold ' ) ; 
ylabel  (  ' Y  (meters)  ' ,  ' font size ' , 18,  ' f ontweight ' ,  ' bold ' )  ; 


figure ; 
hold  all; 

title (' Attempt  to  Trick  Kalman’) 
plot (xm ( 1 ,  : ) ,xm(3,  : )  ,  ' k '  )  ; 
plot (xom ( 1 ,  : ) , xom ( 3 ,  : ) ,  ' b ' ) ; 
plot (xkm ( 1 ,  : ) , xkm ( 3 ,  : ) ,  ' r ' ) ; 

plot (x3m ( 1 , 1 : size (x3m, 2 ) / 2) , x3m ( 3 f 1 : size (x3mf  2 ) /2 ) ,  'g')  ; 

plot (x4am ( 1 , 1 : size (x4am, 2) ) , x4am ( 3 , 1 : size (x4am, 2 ) ) f ' g ' ) ; 

plot (x3m ( 1 f  size (x3mf  2 ) /2  +  l:size (x3mf  2 ) ) ,x3m(3, size (x3mf  2 ) /2  +  l:size (x3mf 

2) ) , 'g') ; 

plot (x4bm ( 1 , 1 : size (x4bm, 2) )  r x4bm (3 , 1 : size (x4bm, 2 ) ) ,  ' g ' ) ; 
legend  ('True  Missile  Track', 'True  Split  Track'); 
set  (gcf ,  'color',  ' w ' )  ; 
grid  on; 

set (  f indob j (gca, ' type ',' line ') ,  'LineWidth',  1.5); 
xlabel ( ' X  (meters )', ' font size ' , 18 , ' f ontweight ' , ' bold ' ) ; 
ylabel  (  ' Y  (meters )',  ' font size ' , 18 ,  ' f ontweight ' ,  ' bold ' )  ; 

figure ; 
hold  all; 

plot (md3 ( 1 : size (md3 , 2 ) / 2 )  ,  ' r '  )  ; 

plot (md3 (size (md3 ,2) /2+l:size (md3 , 2 ) ) , ' b ' ) ; 

set (gcf, 'color', ' w ' ) ; 

legend  ('First  Split  ,  Second  Split'); 
grid  on; 

set (  f indob j (gca,  type ',' line ') ,  'LineWidth',  1.5); 
xlabel ( ' Sample  ( k) ' , ' font size ' , 18 , ' f ontweight ' , ' bold ' ) ; 
ylabel ('Mean  Distance  Error 

(meters )  ' ,  ' font size ' , 18 ,  ' f ontweight ' ,  ' bold ' ) ; 

figure ; 
hold  all; 

plot ( [zeros (1, 100) , md4a] , ' r ' ) ; 

plot ( [zeros (1, 800) , md4b (n(2)-n(l)+l:size (md4a, 2 ) ) ] , ' b ' ) ; 

plot (meml , ' k ' ) ; 

set (gcf, 'color', ' w ' ) ; 
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legend  ('First  Split ',' Second  Split ',' Measurement  Error') 
grid  on; 

set (  f indob j (gca, ' type ',' line ') ,  'LineWidth',  1.5); 
xlabel ( ' Sample  ( k)  ' ,  ' font size '  ,  18 ,  ' f ontweight ' ,  ' bold ' ) ; 
ylabel('Mean  Distance  Error 

(meters ) ' ,  ' font size ' , 18, ' f ontweight ' , ' bold ' ) ; 

figure ; 
hold  all; 
plot (mdo, ' m ' ) ; 
plot (meml , ' k ' ) ; 

set  (gcf , 'color', 'w') ; 

legend  ('Combined  Estimates ',' Measurement  Error'); 
grid  on; 

set (  findobj (gca, ' type ',' line ') ,  'LineWidth',  1.5); 
xlabel ( ' Sample  ( k) ' , ' font size ' , 18 , ' f ontweight ' , ' bold ' ) ; 
ylabel('Mean  Distance  Error 

(meters ) ' , ' font size ' , 18 , ' f ontweight ' , ' bold ' ) ; 

set (  findobj (gca, ' type ',' line ') ,  'LineWidth',  1.5); 


B.  MODEL  ONE  FUNCTION  FILE 

function  [  xlh,Pl,dl,R]  =  modelone (  xlh, x, PI , H, R, znu) 

K1  =  PI *H ' ^ inv ( H*  PI *H '  +  R)  ; 
xlh  =  xlh  +  Kl* (znu  -  H^xlh) ; 

dl^normfH^ (xlh-x) )  ; 

%Covariance  Update 
n=max (size (xlh) ) ; 

Kl 1  =  (eye (n)  -  K1*H) ; 

PI  =  K11*P1*K11'  +  Kl *R*K1 ' ; 

end 

C.  MODEL  TWO  FUNCTION  FILE 

function  [x2h,P2,d2]  =  modeltwo (x2h, x, P2 , znu, H, R) 

K2  =  P2  *H ' * inv (H^P2  '  +  R)  ; 

x2h  =  x2h  +  K2*(znu  -  H^x2h) ; 

d2=norm(H^ (x2h-x) ) ; 

%Covariance  Update 
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end 


n=max (size (x2h) ) ; 

K22  =  (eye (n)  -  K2*H) ; 

P2  =  K22*P2*K22 '  +  K2*R*K2'; 


D.  MODEL  THREE  FUNCTION  FILE 

function  [ PI , P3 , Ps , Ks , zns , xlh, x3h, dl , d3 , itrick]  = 
model three (delta,x,x3,xsh,Q3,Ps,H,R,F3,Ks,Pl,se,znu) 

sigmax=se ; 
sigmay=se ; 
sigmaz=se ; 

zstrue  =  H*x3; 

t=randn ( 3 , 1 ) . * [ sigmax;  sigmay;  sigmaz]/ 
zns  =  zstrue  +  t; 
esm  =  sqrt(t’*t); 

%Generate  stacked  vectors  and  matrices 
Hs  =  [H  zeros (3,6) ; 

zeros  ( 3 , 6 )  H] ; 

Rs  =  [R  zeros  ( 3 , 3 )  ; zeros  ( 3 , 3 )  R] ; 

Fs  =  [F3  zeros (6,6); 

zeros  (6,6)  F3] ; 

Qs  =  [Q3  zeros  (6,6); 

zeros (6,6)  Q3 ] ; 


%Model  3  Target  assignment 
z 1= [ znu; zns ] ; 
z2= [ zns ; znu ] ; 

zlt  =  zl  -  Hs^xsh; 
z2t  =  z2  -  Hs^xsh; 

S3  =  Hs^Ps^Hs 1  +  Rs; 

nl=(zlt) ' ^inv (S3 ) *zlt ; 
n2=(z2t) ' ^inv (S3 ) *z2t ; 

omegal  =  (exp (- (nl/2 ) ) ) / (det (2*pi*S3 ) A ( 1/2 ) ) ; 
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omega2  =  (exp (- (n2/2 ) ) ) / (det (2*pi*S3 ) A ( 1/2 ) ) ; 

%Use  omega  values  to  associate  target  with  measurement 
if  omegal>=omega2 
z  s = z  1  ; 
i  t  r  i  c  k=  0  ; 

else 

zs=z2 ; 
i  t  r  i  c  k=  1 ; 

end 

%Update 

Ks  =  Ps*Hs ’ *inv (Hs*Ps*Hs 1  +  Rs) ; 
xsh  =  xsh  +  Ks* (zs  -  Hs*xsh) ; 

%Convariance  Update 
m=max (size (xsh) ) / 

Kss  =  (eye (m)  -  Ks*Hs) ; 

Ps  =  Kss*Ps*Kss'  +  Ks*Rs*Ks’; 

xlh=xsh (1:6) ; 
x3h=xsh ( 7 : 12) ; 

Pi  =  Ps(l:6,l:6) / 

P3  =  Ps  (7 : 12,7 : 12)  ; 

dl  =  norm (H* (xlh-x) ) / 
d3  =  norm (H* (x3h-x3) ) / 

xsh  =  Fs*xsh; 

Ps  =  Fs*Ps*Fs'  +  Qs; 


end 

E.  MODEL  FOUR  FUNCTION  FILE  WITH  COMPLETE  JDPA  FILTER 
ANALYSIS  OF  ALL  AVAILABLE  MEASUREMENTS 

function  [  x4h, P4 , znu, znc, d4a, d4b, itrick  ]  =  modelfour( 
delta, x4h,x4,P4,Q4,xlh,H,R,se, znc, znu, PI, P3  ) 

%Target  state  transition  matrix 


F  4  = 

[1 

delta 

0 

0 

0 

0; 

0 

1 

0 

0 

0 

0; 

0 

0 

1 

delta 

0 

0; 

0 

0 

0 

1 

0 

0; 

0 

0 

0 

0 

1 

0; 

0 

0 

0 

0 

0 

1]  ; 

%Measurement 

matrix 

H=  [1 

0  0 

0  0 

0; 

0 

0  1 

0  0 

0; 
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0  0  0  0  1  0]; 


sigmax=se ; 
sigmay=se ; 
sigmaz=se ; 
itrick=0 ; 
x4a=x4 (1:6) ; 
zstrue  =  H*x4a; 

t=randn ( 3 , 1 ) . * [ sigmax;  sigmay;  sigmaz]; 
zna  =  zstrue  +  t; 

x4b=x4 (7:12) ; 
zstrue  =  H*x4b; 

v=randn ( 3 , 1 ). * [ sigmax;  sigmay;  sigmaz]; 
znb  =  zstrue  +  v; 

znc= [ zna; znb] ; 

%  JPDA 

znk= [ znu; znc] ; 


if  znk (7  f  1 ) > 0 

zl=znk (1:3) ; 
z2=znk(4:6) ; 
z3=znk(7:9) ; 

xl=xlh; 
x2=x4 (1:6) ; 
x3=x4 (7 : 12) ; 


z  1 1 

=  zl  -  H * x 1 ; 

zl2 

=  zl  -  H*x2; 

zl3 

=  zl  -  H*x3; 

z2 1 

=  z2  -  H * x 1 ; 

z22 

=  z2  -  H*x2; 

z2  3 

=  z2  -  H*x3; 

z31 

=  z3  -  H * x 1 ; 

z32 

=  z3  -  H*x2; 

z33 

=  z 3  -  H*x3; 

SI  = 

H*P1*H'  +  R; 

S2  = 

H*P3*H ’  +  R; 

S3  = 

H*P3*H ’  +  R; 

gii  = 

(exp (-  ( zll )  ’ *inv 

gl2  = 

(exp (- ( zl2 ) ’ *inv 

gl3  = 

(exp (- ( z 1 3 ) ’ *inv 

%Two  split  scenario 


(SI) *zll/2) ) / (det (2*pi*Sl ) 
(SI) *zl2/2) ) / (det (2*pi*Sl) 
(SI) *zl3/2) ) / (det (2*pi*Sl) 


(1/2) ) 
(1/2)  ) 
(1/2) ) 
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g2 1  =  (exp (- (z21) ' *inv(S2) *zll/2) ) / (det (2*pi*S2) A  (1/2) )  ; 
g22  =  (exp (  -  (z22)  ’ *inv(S2) *z22/2) ) / (det (2*pi*S2) A  (1/2) )  ; 
g23  =  (exp (  -  (z23)  ’ *inv(S2) *z23/2) ) / (det (2*pi*S2) A (1/2) ) ; 
g31  =  (exp (- (z31) ’ *inv(S3) *z31/2) ) / (det (2*pi*S3) A (1/2) ) ; 
g32  =  (exp (  -  (z32)  ’ *inv(S3) *z32/2) ) / (det (2*pi*S3) A (1/2) )  ; 
g33  =  (exp (  -  (z33)  ’ *inv(S3) *z33/2) ) / (det (2*pi*S3) A (1/2) ) ; 

gnorm  =  gll  +  gl2  +  gl3  +  g21  +  g22  +  g23  +  g31  +  g32  +  g33 


hi  =  (gll*g22*g33) /gnorm; 
h2  =  (gll*g23*g32 ) /gnorm; 
h3  =  (gl3*g22 *g31 ) /gnorm; 
h4  =  (gl2 *g2 1 *g33 ) /gnorm; 
h5  =  (g31*gl2 *g23 ) /gnorm; 
h6  =  (gl2*g23*g31 ) /gnorm; 

pll  =  hl+h2 ; 
pl2  =  h4+h6 ; 
pl3  =  h3+h5; 
p21  =  h4+h5 ; 
p22  =  hl+h3; 
p23  =  h2+h6; 
p31  =  h3+h6; 
p32  =  h2+h5; 
p33  =  hl+h4 ; 

if  pll>pl2  &&  pll>pl3 
znu=znk (1:3) ; 
itrick=0 ; 

end 

if  pl2>pll  &&  pl2>pl3 
znu=znk (4:6) ; 
itrick=l ; 

end 

if  pl3>pll  &&  pl3>pl2 
znu=znk (7:9) ; 
itrick=l ; 

end 

if  p21>p22  &&  p21>p23 
znc(l:3)=znk(l:3) ; 
itrick=l ; 

end 

if  p22>pll  &&  p22>pl3 
znc (1:3) =znk ( 4 : 6) ; 
itrick=0 ; 

end 

if  p23>pll  &&  p23>pl2 
znc(l:3)=znk(7:9) ; 
itrick=l; 
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end 


if  p31>pl2  &&  p31>pl3 
znc(4:6)=znk(l:3) ; 
i  t  r  i  c  k=  1 ; 

end 

if  p32>pll  &&  p32>pl3 
znc(4:6)=znk(4:6) ; 
i  t  r  i  c  k=  1 ; 

end 

if  p33>pll  &&  p33>pl2 
znc(4:6)=znk(7:9) ; 
i  t  r  i  c  k=  0  ; 

end 

elseif  znk(4,l)>0  %One  split  scenario 


zl=znk (1:3) 

; 

z2=znk (4:6) 

' 

xl=xlh; 

x2=x4 (1:6) ; 

zll  =  z 1  - 

H*xl  ; 

zl2  =  zl  - 

H*x2  ; 

z2 1  =  z2  - 

H*xl  ; 

z2 2  =  z2  - 

H*x2  ; 

SI  =  H^Pl^H 

'  +  R 

S2  =  H*P3*H 

'  +  R 

gll  =  (exp (- (zll)  ’ *inv(Sl) *zll/2) ) / (det (2*pi*Sl) A (1/2) ) 
gl2  =  (exp (- (zl2)  ’ *inv(Sl) *zl2/2) ) / (det (2*pi*Sl) A (1/2) ) 
g2 1  =  (exp (- (z21) ’ *inv(S2) *z21/2) ) / (det (2*pi*S2) A (1/2) ) 
g22  =  (exp (  -  (z22)  ’ *inv (S2) *z22/2) ) / (det (2*pi*S2) A (1/2) ) 

gnorm  =  gll  +  gl2  +  g21  +  g22; 

hi  =  (gll*g22 ) /gnorm; 
h2  =  (g2 l*gl2 ) /gnorm; 

p 1 1  =  hi; 
pl2  =  h2 ; 
p2 1  =  h2 ; 
p22  =  hi; 

if  pll>pl2 

znu=znk  (1:3)  ; 
znc (1:3) =znk(4 : 6) ; 


78 


znc  (4 : 6) =zeros (3, 1)  ; 
itrick=0 ; 

else 

znc ( 1 : 3 ) =znk ( 1 : 3 ) ; 
znu=znk (4:6) ; 
znc  (4 : 6) =zeros (3, 1)  / 
itrick=l ; 

end 

else 

znu=znk (1:3) ; 
znc=znk (4:9) ; 
itrick=0; 

end 

zna=znc (1:3) ; 
znb=znc (4:6) / 

if  x4  (7)  >  0 

x4a=x4  (1:6)  ; 
x4ha=x4h (1:6) / 

P4a=P4 (1:6, 1:6) ; 

%Drop  split  target  estimates  if  >  50  km  from  missile  estimates 
if  norm  (H*  (xlh-x4ha)  )  <  80*10A3 

%  zstrue  =  H*x4a; 

%  t=randn ( 3 , 1 ) . * [ sigmax;  sigmay;  sigmaz]/ 

%  zna  =  zstrue  +  t; 

%Update 

K4a  =  P4a*H ' ^inv (H^P4a^H ’  +  R)  / 
x4ha  =  x4ha  +  K4a* (zna  -  H^x4ha) / 

d4a=norm (H* (x4ha-x4 (1: 6) ) ) ; 

%Covariance  Update 
n=max (size (x4ha) ) / 

K44a  =  (eye (n)  -  K4a*H) / 

P4a  =  K44a^P4a^K44a'  +  K4a*R*K4a1/ 

%Prediction 
x4ha=F4  ^x4ha / 

P4a  =  F4  *P4a*F4 ’  +  Q4; 

else 

P4a=P4a; 

x4ha=x4ha; 

d4a=norm (H* (x4ha-x4 (7:12) ) ) ; 
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end 


x4b=x4 (7:12) ; 
x4hb=x4h (7 : 12) ; 

P4b=P4 (7:12,1:6) ; 

if  norm (H* (xlh-x4hb) )  <  80*10A3 

zstrue  =  H*x4b; 

v=randn ( 3 , 1 ) . * [ sigmax;  sigmay;  sigmaz]; 
znb  =  zstrue  +  v; 

%Update 

K4b  =  P4b*H' *inv (H*P4b*H’  +  R) ; 
x4hb  =  x4hb  +  K4b* (znb  -  H*x4hb) / 

d4b=norm (H* (x4hb-x4 (7:12) ) ) ; 

%Covariance  Update 
n=max (size (x4hb) ) ; 

K44b  =  (eye (n)  -  K4b*H) / 

P4b  =  K44b*P4b*K44b'  +  K4b*R*K4b'; 

%Prediction 
x4hb=F4  *x4hb / 

P4b  =  F  4  *P4b*F4 '  +  Q4; 

else 

P4b=P4b ; 
x4hb=x4hb; 

d4b=norm (H* (x4hb-x4 (7:12) ) ) ; 

end 

elseif  x4 (1)  >  0 

x4a=x4 (1:6) ; 
x4ha=x4h (1 : 6) ; 

P4a=P4 (1:6, 1:6) ; 

P4b=zeros (6,6); 
x4hb=zeros (6,1) ; 
d4b=0 ; 

%Drop  split  target  estimates  if  >  50  km  from  missile  estimates 
if  norm(H* (xlh-x4ha) )  <  80*10A3 

zstrue  =  H*x4a; 

t=randn ( 3 , 1 ). * [ sigmax;  sigmay;  sigmaz]; 
zna  =  zstrue  +  t; 
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%Update 

K4a  =  P4a*H ' *inv (H*P4a*H '  +  R) ; 
x4ha  =  x4ha  +  K4a*  (zna  -  H*x4ha)  ; 

d4a=norm (H* (x4ha-x4 (1: 6) ) ) ; 

%Covariance  Update 
n=max (size (x4ha) ) ; 

K44a  =  (eye (n)  -  K4a*H) ; 

P4a  =  K44a*P4a*K44a'  +  K4a*R*K4a'; 


%Prediction 
x4ha=F4  *x4ha ; 

P4a  =  F4*P4a*F4 '  +  Q4; 

else 

P4a=P4a; 

x4ha=x4ha; 

d4a=norm (H* (x4ha-x4 (1:6)))/ 

end 


else 

x4a=zeros  (6,1)  / 
x4ha=zeros  (6,1)  ; 
zna=zeros  (3,1) / 
d4a=0 ; 

P4a=10 A4 . *  [1 

0 
0 
0 
0 
0 


x4b=zeros  ( 6 
x4hb=zeros  ( 
znb=zeros  ( 3 
P4b=10 A4  .  * 


,D  ; 

6,1)  / 
,D  ; 
[l 
0 
0 
0 
0 
0 


0 

1 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 


0 

0 

1 

0 

0 

0 


0 

0 

1 

0 

0 

0 


0 

0 

0 

1 

0 

0 


0 

0 

0 

1 

0 

0 


0  0; 

0  0; 

0  0; 

0  0; 

1  0; 

0  1]; 


0  0; 

0  0; 

0  0; 

0  0; 

1  0; 

0  1]; 


d4b=0 ; 

end 


%Covariance 


%Covariance 


znc= [ zna; znb] ; 
x4h= [x4ha; x4hb] ; 
P4= [P4a; P4b] ; 


end 


vector 


vector 
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MODEL  PROBABILITY  FUNCTION  FILE  WITH  IMM  MIXING 


function 

[xOlh,  x02h,  P01 ,  P02 ,  ul , u2 , xo, Po, do, clb, c2b] =modelprob (xlh,x2h,Pl,P2,ul,u 
2, znu,H,R,x,clb,c2b) 

%Score  the  Association 
zlt  =  znu  -  H*xlh; 
z2t  =  znu  -  H*x2h; 

51  =  H*P1*H'  +  R; 

52  =  H*P2  *H '  +  R; 

if  sqrt (det (SI ) ) <10 A-5 
scorel  =  10A-5; 

else 

scorel  =  (exp (- (zlt) ’ *inv (SI) *zlt/2 ) ) / (det (2*pi*Sl) A (1/2 ) ) ; 
end 

if  sqrt (det (S2 ) ) <10 A-5 
score2  =  10A-5; 

else 

score2  =  (exp  (- (z2t)  ’ *inv ( S 2 ) *z2t/2) ) / (det (2*pi*S2) A (1/2) )  ; 

end 

%Update  Model  1  and  2  Probabilities 
c  =  scorel*clb  +  score2*c2b; 
ul  =  scorel* (clb/c) ; 
u2  =  score2* (c2b/c) / 

%Combined  Estimates 
xo  =  ul*xlh  +  u2*x2h; 

Po  =  ul* (P1+ [xlh-xo] * [xlh-xo] 1 )  +  u2* (P2+ [x2h-xo] * [x2h-xo] ’ ) ; 

do=norm(H* (xo-x) ) ; 

%Probability  of  changing  state 

pll  =  0.75; 

pl2  =  0.25; 

p2 1  =  0.45; 

p22  =  0.55; 

%Model  1  and  Model  2  mixing 
clb  =  pll*ul  +  p21*ru2; 
c2b  =  pl2^ul  +  p22*u2; 

xOlh  =  xlh* ( (pll*ul) /clb)  +  x2h* ( (p21*u2) /clb) ; 
x02h  =  xlh* ( (pl2*ul) /c2b)  +  x2h* ( (p22*u2 ) /c2b) ; 

ul 1— (pll*ul) /clb; 
u21= (p21*u2 ) /clb; 
u!2= (pl2*ul) / c2b ; 
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u22=(p22*u2) /c2b; 

xlltil=xlh-x01h; 

x21til=x2h-x01h; 

xl2til=xlh-x02h; 

x22til=x2h-x02h; 

P01=ull* [Pl+xlltil*xlltilf ] +u2 1 * [ P2+x2 ltil*x2 1 til ' ]  / 

P02=ul2  * [Pl+xl2til*xl2till ] tu22* [ P2+x22til*x22 til ' ] / 

PO 1=P1 ; 

P02=P2 ; 

end 

G.  PREDICTION  FUNCTION  FILE 

function 

[xlh,  x2h,  PI ,  P2]  ^prediction  (xOlh,  x02h,  P01 ,  P02 ,  FI ,  F2 ,  Q1 ,  Q2 ,  ul ,  u2 ) 

%Predicted  target  motion  after  mixing  of  Model  1  and  2 
xlh  =  Fl*x01h; 
x2h  =  F2*x02h; 

PI  =  F1*P01*F1 ’  +  Ql; 

P2  =  F2*P02*F2'  +  Q2 ; 


end 

H.  STANDARD  KALMAN  FILTER  FUNCTION  FILE  WITH  A  NEAREST 
NEIGHBOR  ALGORITHM  FOR  MEASUREMENT  ASSOCIATION 

function  [  xk,  Pk, dk, ktrick]  =  kalman (  delta , xk, x, PkA znk, Qk, R) 


%Target  state 
FI  =  [1  delta 
0  1 
0  0 
0  0 
0  0 
0  0 


transition  matrix 
0  0  0  0; 

0  0  0  0; 

1  delta  00; 

0  1  0  0; 

0  0  10; 

0  0  01]; 


%Measurement  matrix 
H= [1  0  0  0  0  0 ; 

0  0  1  0  0  0; 

0  0  0  0  1  01  ; 


%Determine  which  measurement  is  closest  to  the  Kalman  estimate 


dl=sqrt ( ( znk ( 1 ) -xk ( 1 ) ) A2  + 
d2=sqrt ( ( znk ( 4 ) -xk ( 1 ) ) ^2  + 
d3=sqrt ( ( znk ( 7 ) -xk ( 1 ) ) ^2  + 


( znk ( 2 ) -xk ( 3 ) ) A2) ; 
( znk ( 5 ) -xk  ( 3 )  ) A2)  ; 
( znk ( 8 ) -xk ( 3 ) ) A2) ; 
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ktrick=0 ; 


if  abs  ( znk ( 4 ) ) >0  &&  abs ( znk ( 7 ) ) >0 

if  d2<dl  &&  d2<d3 
znc=znk (4:6) ; 
ktrick=l ; 

elseif  d3<dl  &&  d3<d2 
znc=znk (7:9) ; 
ktrick=l ; 
else 

znc=znk(l:3) ; 
ktrick=0 ; 

end 

elseif  abs ( znk ( 4 ) ) >0 
if  d2<dl 

znc=znk (4:6) / 
ktrick=l ; 
else 

znc=znk (1:3) / 
ktrick=0 ; 
end 


else 

znc=znk ( 1 : 3 ) ; 
ktrick=0 ; 

end 

%Update 

Kk  =  Pk*H ' ^inv (H^Pk^H '  +  R) ; 
xk  =  xk  +  Kk* (znc  -  H*xk) ; 

dk  =  norm (H* (xk-x) ) ; 

%Covariance  Update 
n=max (size (xk) ) / 

Kkk  =  (eye (n)  -  Kk*H) ; 

Pk  =  Kkk*Pk*Kkk '  +  Kk*R*Kk ’ ; 

%Predicted  Target  Motion 
xk  =  Fl*xk; 

Pk  =  Fl*Pk*Fl 1  +  Qk; 

end 
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